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Abstract 

We  develop  a  9.6  million  node  unstructured  grid  finite-element  forward  barotropic  model  in  order 
to  understand  the  shallow  water  dynamics  and  dissipation  mechanisms  of  the  Indian  and  Western 
Pacific  Oceans  down  to  sub-kilometer  scale  at  the  coast.  Tidal  model-data  comparisons  against  tide 
gauges  and  a  global  data-assimilated  tidal  model  are  conducted  in  order  to  identify  the  capabilities 
and  limitations  of  our  model.  The  average  root-mean-square  (RMS)  discrepancies  of  the  total  free 
surface  at  coastal  tide  gauges  is  14  cm,  ~3  cm  smaller  than  the  data-assimilated  model.  Sensitivities 
related  to  lateral  boundary  conditions,  bathymetry,  and  dissipative  processes  are  explored.  Lateral 
boundary  conditions  were  found  to  induce  global  resonant  effects  on  the  lunar  semi-diurnal  modes 
when  poorly  placed  elevation  specified  boundary  conditions  were  used.  This  problem  is  largely 
resolved  by  using  an  absorption-generation  layer  at  the  boundary,  although  errors  may  persist  in 
the  flux  component  obtained  from  data-assimilated  tidal  models.  Parameterization  of  internal  tide 
dissipation  is  identified  as  the  most  important  aspect  to  control  deep  water  solutions,  and  reduce 
the  RMS  discrepancies  of  the  entire  system.  Two  forms  of  this  parameterization  were  presented  and 
their  spatial  distributions  of  dissipation  were  compared.  Bathymetry  has  a  negligible  effect  on  the 
tidal  solutions  in  deep  water,  but  local  high-resolution  bathymetry  results  in  significant  reductions  to 
the  average  RMS  discrepancies  on  the  continental  shelf  (26%)  and  at  the  coast  (30%).  Implementing 
a  spatially  varying  bottom  friction  coefficient  based  on  sediment  types  decreases  the  average  RMS 
discrepancy  at  the  coast  by  9%  predominantly  due  to  its  positive  effects  in  the  Yellow  Sea.  The 
model  is  shown  to  capture  a  large  amount  of  the  tidal  physics  and  has  the  potential  for  application 
to  a  wide  range  of  shallow  water  problems  in  the  region. 

Keywords:  finite-element;  unstructured  grid;  barotropic  tides;  bathymetry;  internal  tide  genera¬ 
tion;  bottom  friction 


1  Introduction 


The  Indian  Ocean  and  Western  Pacific  Ocean  (IndWPac)  represent  approximately  30%  of  the  surface 
area  of  the  world  ocean,  are  interconnected  by  marginal  seas  such  as  the  Java,  Timor,  Banda, 
Andaman  and  Arafura  Seas,  and  are  separated  by  the  intricate  island  chains  of  Indonesia  and  the 
Philippines.  Major  ports  and  cities  are  located  in  the  northern  parts  of  both  the  Indian  Ocean 
(Dubai,  Karachi,  Mumbai,  Colombo)  and  the  Western  Pacific  Ocean  (Hong  Kong,  Shanghai,  Tokyo, 
Singapore),  representing  a  significant  portion  of  the  world’s  economy  and  human  population. 

Given  the  rapidly  changing  coastal  zone  in  much  of  the  IndWPac  region,  the  development  of 
a  large-scale  comprehensive  barotropic  model  is  needed  to  advance  modeling  capabilities  of  tides, 
wind  driven  surge  and  circulation  and  their  nonlinear  interactions  nearshore.  We  aim  to  be  on 
par  with  contemporary  models  in  other  regions  receiving  significant  attention  such  as  the  western 
North  Atlantic  (e.g.  Westerink  et  al.,  2008;  Bunya  et  al.,  2010;  Kerr  et  al.,  2013).  The  challenge  of 
the  IndWPac  region  is  that  it  is  more  complex  in  its  geometry,  topography  (e.g.  it  contains  many 
interconnected  shallow  seas  and  island  chains)  and  associated  hydrodynamics  than  the  Western 
North  Atlantic  region. 

To  model  the  dynamics  at  coastal  locations  within  the  IndWPac  region,  all  processes  and  ex¬ 
changes  from  ocean  basin  scale  to  harbor  inlet  scale,  must  be  appropriately  represented.  Coarse 
resolution  global  models  (e.g.  Egbert,  Ray,  and  Bills,  2004;  Green  and  Nycander,  2013;  Buijsman 
et  ah,  2015)  capture  the  large-scale  dynamics  but  they  may  miss  nearshore  features  and  important 
nonlinearities.  Conversely,  carefully  calibrated  local  domain  hydrodynamic  models  can  accurately 
capture  local  effects  (e.g.  Green  and  David,  2013;  Zu,  Gan,  and  Erofeeva,  2008,  in  South  China  Sea), 
but  they  may  miss  basin- scale  effects  from  offshore  and  are  restricted  by  boundary  conditions.  A 
single  unstructured  grid  finite-element  model  of  the  region  is  suitable  because  it  can  utilize  varying 
resolutions  to  produce  high  fidelity  coastal  bathymetry  of  critical  geographic  and  topographic  fea¬ 
tures  such  as  island  chains,  reef  systems,  and  floodplain  systems;  provide  connectivity  to  estuarine 
and  harbor  systems  where  most  water  level  gauges  are  located  and  where  dense  coastal  populations 
live;  and  capture  key  dynamics  of  a  large  regional  domain. 

To  assess  the  capabilities  and  limitations  of  a  large-scale  finite-element  barotropic  model  in 
the  IndWPac  region,  model-data  comparisons  of  the  tidal  elevations  are  investigated  in  this  study. 
This  can  be  achieved  point-wise  at  tide  gauges  or  regionally  against  global  data-assimilative  model 
atlases.  Examples  of  the  latter  include  TPXO  (Egbert  and  Erofeeva,  2002),  FES  (Lyard  et  al., 
2006),  and  NA0.99b  (Matsumoto,  Takanezawa,  and  Ooe,  2000).  These  models  assimilate  data 
from  satellite  altimetry  and  selected  coastal  tide  gauges  to  accurately  obtain  estimates  of  the  tidal 
elevation  fields  in  terms  of  individual  harmonic  constituents.  Modern  data-assimilated  models  can 
achieve  M2  tidal  wave  root-mean-square  errors  (RMSE)  of  0.5  -  0.7  cm  versus  deep-ocean  bottom 
pressure  recorder  stations  (Stammer  et  al.,  2014).  In  contrast,  M2  RMSE  ranges  between  5.6  - 
12.7  cm  for  purely  hydrodynamic  global  models  without  data-assimilation  (Stammer  et  al.,  2014). 
However,  non-assimilative  forward  models  on  large-domains  can  be  applied  to  a  wide  variety  of 
problems  including  wind,  pressure  and  wave  coupling  effects,  and  may  be  used  to  conduct  future 
forecasting  and  perturbation  response  analysis  (Green  and  David,  2013),  e.g.  due  to  changing  sea 
level,  dredging  operations,  and  land  reclamation  (Suh,  Lee,  and  Kim,  2014). 

Based  on  the  model-data  comparisons  we  aim  to  discern  areas  where  a  high-resolution  model  has 
advantages  over  coarser  models  and  where  lingering  discrepancies  remain.  Moreover,  a  discussion 
on  the  major  reasons  for  these  is  of  great  interest  to  inform  future  developments  of  coastal  ocean 
models  and  their  applications.  This  is  achieved  through  assessment  of  the  sensitivities  of  the  various 
controls  on  the  barotropic  tidal  dynamics  (lateral  boundary  position,  boundary  condition  imple¬ 
mentation,  bathymetry,  and  dissipative  processes).  For  example,  two  global  bathymetric  databases 
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are  compared,  and  high- resolution  local  bathymetric  data  is  included  where  available  to  assess  the 
potential  of  improved  bathymetry  facilitating  improvements  in  the  solution. 

Regarding  dissipative  processes,  significant  progress  has  been  made  in  reducing  tidal  elevation 
errors  in  deep  water  through  parameterizations  of  tidal  energy  conversion  from  barotropic  to  baro- 
clinic  modes  through  internal  tide  generation  over  submarine  ridges  (Jayne  and  St.  Laurent,  2001; 
Egbert,  Ray,  and  Bills,  2004;  Zaron  and  Egbert,  2006;  Green  and  Nycander,  2013;  Buijsman  et  al., 
2015).  This  paper  will  refer  to  that  physical  process  as  internal  tide  dissipation ,  and  will  discuss 
the  significance  of  its  effects  in  the  IndWPac  region.  In  addition,  spatially  varying  bottom  friction 
coefficients  are  rarely  considered  in  large-scale  models.  Instead,  a  canonical  spatially  constant  coef¬ 
ficient  is  commonly  applied  (Lyard  ct  al.,  2006;  Egbert  and  Erofeeva,  2002).  However,  changing  the 
bottom  friction  coefficient  has  been  shown  to  have  positive  effects  regionally  (Lefevre,  Provost,  and 
Lyard,  2000).  We  take  a  preliminary  look  at  the  effect  of  spatially  varying  coefficients  that  are  based 
on  local  sediment  types  and  hydrodynamics,  an  important  consideration  for  the  IndWPac  domain 
where  diverse  sediment  categories  are  present  across  wide  and  very  shallow  continental  shelves. 

To  summarize,  this  paper  describes  the  development  of  the  IndWPac  unstructured  grid  and  hy¬ 
drodynamic  modeling  system,  built  with  state-of-the-art  bathymetric  datasets,  absorptive  boundary 
conditions,  and  parameterizations  of  internal  tide  and  bottom  friction  dissipation.  We  analyze  the 
sensitivity  of  the  model  to  these  four  factors,  and  conduct  model-data  comparisons  against  both 
tide  gauge  records  and  a  data-assimilated  tidal  model.  The  capabilities  and  limitations  of  the  model 
are  identified  and  discussed.  Suggested  areas  of  focus  in  order  to  advance  barotropic  coastal  ocean 
models  are  highlighted. 


2  Domain  Definition,  Bathymetry,  and  Unstructured  Grid 
Development 


Our  basin-scale  model  includes  the  entire  Indian  Ocean,  the  western  half  of  the  Pacific  Ocean,  and 
the  Southern  Ocean  between  these  extents.  Specifically,  the  domain  (Fig  1)  lies  between  17.9°E  - 
175. 8°E  longitude  and  73.3°S  -  62.7°N  latitude  covering  an  area  of  roughly  150  million  km2.  There 
are  two  open  ocean  boundaries:  a  longitudinal  parallel  boundary  running  from  nearby  the  Cape 
of  Good  Hope,  South  Africa  to  Antarctica;  and  a  concave  shaped  boundary  between  the  Bering 
Sea  coast  of  Kamchatka  Krai,  Russia  and  Antarctica.  The  boundaries  were  chosen  so  that  tidal 
amphidromic  points  and  complications  with  the  Aleutian,  Hawaiian  and  New  Zealand  islands  in  the 
Pacific  Ocean  were  avoided  (an  illustration  on  the  effects  of  boundary  placement  is  shown  in  §5.1). 

The  mesh  is  a  triangular  unstructured  grid  with  resolution  ranging  from  as  large  as  25  km  in 
parts  of  the  deep  ocean  down  to  1  km  along  most  coastlines  (Fig  2).  Additionally,  resolution  is  as 
fine  as  100  m  in  the  ports  and  harbors  of  Hong  Kong,  Tokyo  Bay  and  Osaka  Bay.  The  grid  contains 
a  total  of  9.6  million  nodes  and  18.8  million  elements. 

Development  of  the  unstructured  grid  was  achieved  predominantly  through  an  automated  al¬ 
gorithm  that  we  developed  in-house  based  on  the  MATLAB  DistMesh  code  (Persson  and  Strang, 
2004).  Resolution  is  varied  through  an  edgelength  (local  grid  resolution)  function  A defined  as  the 
minimum  of  three  criteria: 


A/.;  =  min 


Am  4-  cktid. 


2  IT  h  \ 
a,  \Vh\) 


(1) 


where  A m  is  the  nominal  minimum  edgelength,  d  is  the  distance  from  a  node  to  the  closest  coast¬ 
line  boundary,  T  is  the  period  of  the  M2  tidal  wave,  h  is  the  bathymetric  depth,  and  a  are  the 
dimensionless  user-defined  coefficients  for  each  criterion:  distance  from  the  coastline  (a^  =  0.075), 
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Figure  1:  Bathymetric  depths  of  the  grid  as  interpolated  from  various  sources  (Table  1)  using  a  cell-averaged 
approach. 


wavelength  (a^  =  600),  and  topographic  length  scale  (a5  =  30).  These  criteria  ensure  that  im¬ 
portant  bathymetric  features  are  captured  throughout  the  ocean  in  addition  to  obtaining  higher 
resolution  nearshore. 

Model  bathymetry  is  interpolated  onto  the  grid  from  a  number  of  sources  in  a  specified  order 
using  an  automated  cell-averaging  technique  (Bilskie  and  Hagen,  2013)  as  summarized  in  Table  1. 
The  adopted  background  bathymetry  is  the  1/120°  SRTM30.PLUS  global  database  combined  with 
a  synthetic  realization  of  seafloor  roughness  along  the  abyssal  hills  (Goff  and  Arbic,  2010;  Timko 
et  al.,  2017).  The  synthetic  abyssal  hill  roughness  is  used  because  the  effective  resolution  of  the 
global  altimetric  based  bathymetry  is  limited  to  >10  km  in  the  deep  ocean  while  ~1  km  resolution 
is  necessary  to  describe  the  required  topographic  roughness  that  generates  internal  tides  converting 
barotropic  energy  into  baroclinic  energy  (Goff  and  Arbic,  2010;  Melet  et  ah,  2013;  Timko  et  al.,  2017). 
In  addition,  to  include  depths  under  ice  shelves  in  Antarctica  we  interpolate  from  the  TPX08  model 
bathymetry  containing  the  Padman  et  al.  (2002)  dataset. 

For  shallower  regions  (in  depths  <  500  m)  where  the  abyssal  hill  roughness  is  not  impor¬ 
tant,  we  start  by  interpolating  from  the  global  1/240°  SRTM15JPLUS  database  which  improves 
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Figure  2:  Resolution  of  the  unstructured  grid,  which  varies  based  on  topographic  gradients,  depths  and 
proximity  to  the  coastline.  The  resolution  at  the  coastline  is  ~1  km  in  most  regions,  and  up  to  ~25  km  in 
the  deep  and  flat  regions  of  the  ocean. 


on  SRTM30-PLUS  with  newer  measured  nearshore  bathymetry  and  topography  sources  thereby  re¬ 
ducing  the  number  of  erroneous  holes  in  the  data.  On  top  of  this,  100  m  Deepreef  Explorer  Great 
Barrier  Reef  and  Coral  Seas  (GBR),  and  Kerguelen  Plateau  (KP)  datasets  were  applied.  It  was 
discovered  that  Deepreef  Explorer  GBR  in  the  Torres  Strait/New  Guinea  Region  matches  substan¬ 
tially  better  with  GEBCO_2014  than  SRTM15-PLUS,  thus  GEBCO_2014  was  applied  locally  here 
(differences  between  the  two  databases  are  discussed  further  in  §5.2).  Also,  90  m  East  Asia  nearshore 
bathymetry  datasets  in  the  Philippines,  Japan,  Gulf  of  Thailand,  South  China  Sea,  and  East  China 
Sea  regions;  and  local  high-resolution  bathymetry  and  grids  privately  obtained  for  Tokyo  Bay  and 
South  Korea  were  applied.  However,  even  in  the  high-resolution  datasets,  erroneous  depth  in  harbor 
complexes  and  channels  persist.  These  were  corrected  where  possible  using  data  from  FUGAWI 
navigational  charts  (https://www.fugawi.com/).  However,  the  errors  in  the  final  bathymetry  that 
is  applied  to  IndWPac  are  still  largely  uncertain. 
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3  ADCIRC  Hydrodynamic  Model 

3.1  Governing  Equations 

The  horizontal  two-dimensional  implementation  of  the  Advanced  Circulation  coastal  ocean  model 
(ADCIRC-2DH)  is  used  to  calculate  the  hydrodynamics  (Westerink  et  al.,  2008;  Westerink  et 
al.,  1992).  The  governing  equations  are  the  shallow  water  equations  (SWE)  in  primitive,  non¬ 
conservative,  and  barotropic  form: 

Yt  +  V  '  +  cr(a!)(T/  -  7lc)  =  o  (2) 

du  _  x  _  lulu  _ 

—  +  u-Vu  +  /kxu  +  ^V(77  -  T]EQ  -  t]sal)  +  <?/— —  +  Cu 

H  (3) 

—  77V  •  [i/tH(Vu  +  VuT)]  +  <j{x)(u  —  uc)  =  0 

where  77  is  the  surface  elevation,  H  =  h  +  77  is  the  total  water  depth  in  which  h  is  the  still  water 
depth,  u  is  the  depth- averaged  velocity  vector,  g  is  the  acceleration  due  to  gravity,  k  is  the  vertical 
unit  vector,  and  f  —  20.  sin  <fi  is  the  Coriolis  parameter  in  which  O  is  the  angular  speed  of  the  earth, 
and  <j>  is  the  latitude.  The  quantity  t}eq  is  the  equilibrium  tide,  and  t)sal  is  the  ocean  self- attraction 
and  loading  term  (SAL).  In  the  dissipation  terms,  Cf  is  the  coefficient  of  bottom  friction,  C  is  the 
dissipation  matrix  due  to  the  generation  of  internal  tides,  and  vt  is  the  horizontal  eddy  viscosity 
coefficient.  Finally,  we  impose  an  absorption-generation  sponge  layer  (e.g.  Zhang  et  al.,  2014)  where, 
<j{x)  are  the  spatially  varying  absorption  coefficients  applied  over  the  defined  sponge  boundary,  and 
7/c  and  uc  are  the  corresponding  reference  solutions  for  surface  elevation  and  velocity  respectively. 

3.2  Ocean  Self-attraction  and  Loading  Term 

The  ocean  self- attracting  and  loading  (SAL)  term,  t]sal  is  related  to  the  yielding  of  the  solid  Earth 
to  tides  and  to  the  weight  of  the  ocean  and  its  self-attraction  (Hendershott,  1972).  For  the  large  scale 


Table  1:  Bathymetric  data  sources,  location  applied,  resolution  and  availability.  Interpolation  onto  our 
grid  is  conducted  in  the  order  shown  in  this  table 


Name 

Source  (s) 

Location 

Resolution 

Availability 

SRTM30.PLUS 

Becker  et  al.  (2009) 

globally  >500  m  depth 

1/120° 

free  at  website1 

Abysall  Hills 

Goff  and  Arbic  (2010)  and  Melet  et  al.  (2013) 

globally  >500  m  depth 

1/120° 

prvt.  comm. 

SRTM15_PLUS 

Sandwell  et  al.  (2014) 

globally  <500  m  depth 

1/240° 

free  at  website2 

TPX08 

Padman  et  al.  (2002) 

<65°S 

1/30° 

free  at  website3 

GEBCO-2014 

Weatherall  et  al.  (2015) 

Torres  Strait/New  Guinea 

1/120° 

free  at  website4 

Deepreef  Explorer  GBR 

Beaman  (2010) 

Great  Barrier  Reef  &  Coral  Sea 

1/1000° 

free  at  website5 

Deepreef  Explorer  KP 

Beaman  and  O’Brien  (2011) 

Kerguelen  Plateau 

1/1000° 

free  at  website6 

TCarta  Marine 

TCarta  Marine  (2012) 

East  Asia  nearshore 

1/1200° 

proprietary7 

Tokyo  Bay  HR 

Shintaro  Bunya  ( prvt .  comm.,  2015) 

Tokyo  Bay 

FE  grid 

prvt.  comm. 

South  Korea  HR 

SeungWon  Suh  (prvt.  comm.,  2017) 

South  Korea 

FE  grid 

prvt.  comm. 

Harbor  hand-edits 

FUGAWI  Navigational  Charts 

various  harbors  and  channels 

FE  grid 

- 

FE  grid:  indicates  data  was  received  on  a  finite-element  grid 
1 :  ftp : //t opex . ucsd . edu/pub/srtm30_plus/ 

2:  ftp : //topex . ucsd.edu/pub/srtml5_plus/ 

3:  http : //volkov.oce .orst . edu/tides/tpxo8_atlas .html 
4 :  http : //www . gebco . net/data_and_products/gridded_bathymetry_data/ 
5:  https : / /www . deepreef . org/bathymetry/65-3dgbr-bathy . html 
6:  https : / /www . deepreef . org/bathymetry/98-kergdem-bathy . html 
7:  provided  by  Factory  Mutual  Insurance  (FM  Global),  Norwood,  MA 
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IndWPac  domain  it  is  essential  to  include  the  effect  of  SAL  terms  on  the  tides.  However,  since  the 
model  is  regional,  the  global  integrals  of  the  tidal  elevations  required  to  be  solved  iteratively  for  the 
SAL  terms  (Ray,  1998)  are  not  available.  Thus,  in  this  study  the  amplitudes  and  phases  of  SAL  for 
each  tidal  constituent  are  simply  interpolated  from  those  used  in  the  global  data-assimilated  model 
FES2014  (Lyard  et  al.,  2006)  onto  our  mesh  and  forced  by  reconstructing  the  time  series  from  the 
constituents.  Given  the  accuracy  of  state-of-the-art  global  data-assimilated  models  (Stammer  et  al., 
2014),  the  slowly  varying  SAL  terms  obtained  from  these  models  are  also  assumed  to  be  sufficiently 
accurate.  However,  the  calculation  of  SAL  through  global  integrals  to  obtain  full  consistency  with 
the  surface  elevation  (including  non-periodic  components)  is  ultimately  desired  (c.f.  Apecechea  et  al., 
2017). 


3.3  Dissipation  due  to  the  Generation  of  Internal  Tides 


Internal  tides  generated  by  flow  over  rough  bathymetry  are  major  contributors  to  barotropic  tidal 
energy  dissipation  in  the  deep  ocean,  equivalent  to  around  30%  of  the  global  total  (Egbert  and  Ray, 
2001).  As  a  result,  parameterization  of  this  dissipation  is  necessary  in  barotropic  numerical  models 
that  include  expanses  of  ocean  where  major  submarine  ridges,  island  chains  and  shelf  breaks  that 
induce  internal  waves  are  present.  In  this  study,  parameterization  of  internal  tide  dissipation  is 
particularly  important  since  the  Indian  Ocean  basin  contains  narrow  shelves  and  vast  expanses  of 
open  ocean  where  the  dissipation  due  to  internal  tides  over  its  well  defined  abyssal  hills  are  crucial 
to  the  accuracy  of  the  tidal  solutions. 

Parameterizations  of  internal  tide  dissipation  are  usually  based  on  a  linear  wave  drag  type  imple- 

1  /2 

mentation  valid  only  for  subcritical  topography  (7  <  1).  Here,  7  =  ,  in  which  a  = 

is  the  internal  wave  slope,  uj  is  the  angular  frequency  of  the  pertinent  tidal  wave  (M2  in  this  study), 
and  Nb  is  the  Brunt-Vaisala  frequency  at  the  seabed.  In  this  study,  we  investigate  two  subcritical 
theory  parameterizations  for  the  dissipation  matrix  C  in  (3):  one  based  only  on  local  topographic 
features,  and  another  that  includes  the  nonlocal  effects  on  wave  generation. 

First,  we  use  a  simple  and  robust  parameterization  that  takes  into  account  the  directionality  of 
dissipation  (which  we  denote  as  the  ‘ Directional 5  method)  similar  to  that  presented  by  Lyard  et  al. 
(2006)  is: 


C  =  Cm 


\{N2-uj2)(N2-uj2)]1'2 


UJ 


h2x 


hX  hy 

hl . 


(4) 


where  C^ir  is  a  scale  factor,  N  is  the  depth-averaged  Brunt-Vaisala  frequency,  and  the  subscripts 
cx5  and  V  indicate  gradients  in  the  longitudinal  and  latitudinal  directions  respectively.  Note  that  we 
have  substituted  the  typical  wavenumber,  k  in  Lyard  et  al.  (2006)  for  the  fundamental  internal  mode 
at  the  pertinent  tidal  frequency  (Zaron  and  Egbert,  2006).  The  Directional  method  only  dissipates 
across  slopes  (rather  than  along  them). 

Second,  a  rigorous  formulation  for  C  that  includes  the  nonlocal  effects  of  the  nearby  topography 
on  internal  tide  generation  (Melet  et  al.,  2013)  was  derived  by  Nycander  (2005)  (denoted  as  the 
‘Nycander'  method  hereafter).  It  has  the  following  form  in  a  general  coordinate  system  (Green  and 
Nycander,  2013): 


c  =  CAryc^V1-^ 


p 

uj^ 


2  JxK 
J  h*  4-  J  h* 

U Xfly  T-  'Jyll'x 


Jxhy  +  Jyhl 

2 Jyhiy 


(5) 


where  Cnvc  is  a  scale  factor,  and  J  is  a  convolution  integral  of  a  filtered  Green’s  function  of  the 
topographic  heights  h*  (defined  positive  from  seabed)  within  a  specified  radius  from  the  point 
of  interest  (c.f  Green  and  Nycander,  2013;  Nycander,  2005).  Readers  should  refer  to  Green  and 
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Nycander  (2013)  for  the  equations  and  numerical  scheme  to  compute  the  gradients  of  J  (Jx  and  Jy) 
on  the  structured  grid  of  the  bathymetric  data. 

The  advantages  of  the  Directional  method  is  that  C  is  positive  definite,  and  it  does  not  require  the 
computationally  intensive  calculation  of  J  allowing  it  to  be  quickly  implemented  into  the  model.  On 
the  other  hand,  the  Nycander  method  allows  for  the  nonlocal  generation  of  internal  tides.  However, 
C  in  (5)  is  not  guaranteed  to  be  positive  definite  since  the  sign  of  the  gradients  of  J  and  h*  do 
not  necessarily  conform.  Furthermore,  the  calculation  of  the  gradients  of  J  are  computationally 
expensive  so  it  is  not  as  readily  implemented  into  a  numerical  model. 

Note  that  for  h  <  100  m  we  set  C  =  0  in  both  methods,  because  the  topographic  gradients  on 
the  continental  shelf  should  be  small,  and  bottom  friction  dissipation  starts  to  dominate  here. 

Modifications  to  get  a  positive  definite  C,  and  to  incorporate  the  nonlocal  effects  of  iV&  into  the 
Nycander  method  are  implemented  and  briefly  evaluated  in  this  study.  We  also  investigate  whether 
the  nonlocal  wave  generation  effects  captured  by  the  Nycander  method  have  a  meaningful  effect  on 
the  solution  by  comparing  the  results  between  the  twro  methods  (see  §5.3). 

3.4  Bottom  Friction  Dissipation 

Dissipation  due  to  bottom  friction  is  known  to  account  for  a  significant  proportion  of  dissipation  of 
the  barotropic  tides,  particularly  in  shallow  regions  (h  100  m).  It  is  commonly  assumed  that  the 
bottom  stress  is  proportional  to  square  of  the  depth-averaged  velocity  in  barotropic  tidal  simulations 
(e.g.  Buijsman  et  al.,  2015;  Green  and  David,  2013;  Lyard  et  al.,  2006)  as  adopted  in  (3).  Values 
for  the  coefficient  of  bottom  friction,  Cf  have  been  shown  to  range  between  orders  10~3  and  10-2 
based  on  measurements  in  continental  shelf  and  estuarine  regions  (e.g.  You,  2005;  Heathershaw, 
1979;  Heathershaw  and  Simpson,  1978;  Charnock,  1959).  Canonical  global  values  of  Cf  equal  to 
2.5xl0-3  (Lyard  et  al.,  2006)  or  3.0xl0-3  (Egbert  and  Erofeeva,  2002)  are  usually  applied  as  a 
spatial  constant  in  large-scale  tidal  models. 

It  has  been  suggested  that  deviations  from  the  canonical  value  of  Cf  globally  do  not  significantly 
change  the  overall  dissipation  but  that  deviations  by  an  order  of  ten  can  significantly  degrade  the 
tidal  solution  (Lyard  et  al.,  2006).  Nevertheless,  if  other  dissipation  mechanisms  are  reliable  (internal 
tide)  there  is  evidence  that  local  variations  in  Cf  over  the  range  of  physically  plausible  values  (10-3 
to  10-2)  can  improve  local  tidal  solutions  (e.g  Lefevre,  Provost,  and  Lyard,  2000).  In  this  study,  we 
present  a  semi-data-informed  method  of  calculating  spatially  varying  Cf.  The  goal  is  not  to  obtain 
the  ‘optimal’  Cf  but  merely  to  show  that  it  is  possible  to  calculate  spatially  varying  Cf  that  locally 
improves  tidal  solutions  based  on  some  knowledge  of  the  seabed  and  physical  properties  of  the  flow. 

We  start  with  the  log-law  formulation  of  Cf  (Schlichting,  1979): 

Cf  =  [k/1ii(0.5F/zo)]2  (6) 

where  k  —  0.4  is  the  von  Karman  constant,  and  zo  is  the  seabed  roughness  length  which  can  be 
equated  to  an  effective  sediment  roughness,  ks  (=  30zo ).  It  is  important  to  note  that  ks  is  not 
simply  a  function  of  the  sediment  roughness  (grain-size)  itself,  rather  it  is  mainly  determined  by  the 
heights  of  ripples  and  dunes  (bedforms)  that  form  due  to  the  prevailing  currents  which  can  be  a 
major  source  of  the  resultant  bed  stress  (Heathershaw,  1979).  To  estimate  ks  that  takes  into  account 
the  bedform  heights,  we  use  empirical  equations  (Rijn,  2007)  that  are  a  function  of  median  sediment 
grain  diameter  cfeo,  sediment  density  relative  to  water  s,  an  effective  mean  current  speed  u/,  and 
the  depth  h  (see  A).  The  empirical  equations  return  small  values  of  ks  when  either  the  sediments 
are  light  and  the  tidal  currents  are  strong  flattening  out  the  bed,  or  when  the  sediment  grains  are 
too  heavy  for  the  currents  to  create  bedforms.  In  between  these  extremes,  ripples  and  dunes  will 
form  resulting  in  larger  values  of  ks. 
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To  obtain  the  sediment  grain  sizes  we  make  use  of  a  database  of  the  census  of  the  world’s  seafloor 
sediment  types  (Dutkiewicz  et  al.,  2015).  We  map  these  sediment  types  onto  physically  reasonable 
values  of  d$ o  (see  Table  2).  For  pelagic  type  sediments  (oozes  and  clays)  Cf  is  set  to  2.5 xlO-3  as 
a  default  roughness.  Relative  sediment  density  s  =  1.722  (dry  bulk  density  by  mass  of  sand,  Rijn, 
2007)  for  d5o  >  dsand ,  s  =  1.2  (natural  sediment  with  organic  materials  involved,  Rijn,  2007)  for 
c/50  <  dsnt ,  and  is  linearly  interpolated  in  between.  Here,  dsanc[  =  6.2xl0~5  m,  and  dsnt  =  3.2xl0-5 
m,  where  the  assumption  is  made  that  the  finer-sized  sediments  in  the  census  database  contain  a 
higher  percentage  of  lighter  organic  material.  The  effective  mean  current  speed  Uf  is  defined  as 
(Zaron,  2017): 

U/=^  +  0.5^|t7fc|2j  (7) 

where  uq  is  a  constant  non-tidal  current  (Snyder,  Sidjabat,  and  Filloux,  1979),  that  we  set  equal  to 
0.25  m/s  (Zaron,  2017),  and  Uk  (both  x  and  y  components)  is  the  amplitude  of  the  kth  tidal  current 
constituent.  The  spatially  constant  Cf  —  2.5  xlO-3  simulation  is  used  to  approximate  Uf  in  order 
to  compute  the  spatially  varying  Cf  map  (Fig.  3). 

3.5  Lateral  Eddy  Viscosity 

To  include  lateral  eddy  viscosity  in  ADCIRC,  the  Smagorinsky  model  (Smagorinsky,  1963)  is  used 
to  evaluate  vt  since  it  changes  with  the  local  mean  strain-rate  providing  a  degree  of  physical  fidelity, 
as  well  as  providing  stability  to  the  advection  terms.  vt  is  calculated  as  follows  (Dresback,  Kolar, 
and  Luettich,  Jr.,  2005): 

Vt  =  C^Ae^/ (V  •  uT)2  +  (V  x  uT)2  (8) 

where  C M  is  a  free  coefficient  that  we  set  equal  to  0.2,  and  Ae  is  the  local  element  area. 

3.6  Lateral  Boundary  Conditions 

Lateral  open  ocean  boundaries  are  forced  by  reconstructing  the  elevations  from  the  tidal  constituents 
obtained  from  a  global  data- assimilative  model,  TPX08  (Egbert  and  Erofeeva,  2002).  In  this  study 
we  force  with  the  major  semi-diurnal  (M2,  N2,  S2,  K2)  and  diurnal  (Ki,  Oi,  Pi,  Qi)  constituents, 
which  are  also  used  to  force  the  SAL  and  equilibrium  potential  terms.  Prescribing  the  elevations 
at  the  open  boundaries  is  a  reflecting  boundary  condition  that  allows  the  velocities  to  freely  satisfy 
the  governing  equations.  In  some  cases  this  condition  can  generate  spurious  modes  that  may  lead 
to  instabilities.  Utilizing  an  absorption-generation  sponge  layer  can  reduce  the  production  of  these 
modes,  as  demonstrated  in  §5.1. 


Table  2:  Median  grain  sizes  c/50  and  relative  density  s  for  each  sediment  type  used  in  the  calculation  of  Cf 


Sediment  Type 

d50  [m] 

s 

Gravel  and  coarser 

3.0xl0-3 

1.722 

Sand 

l.OxlCT4 

1.722 

Silt 

5.0xl0~5 

1.513 

Ash  and  volcanic  sand/gravel 

l.OxlCT3 

1.722 

Siliceous  mud 

4.0xl0-5 

1.339 

Fine-grained  calcareous  sediment 

4.5xl0-5 

1.426 
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Figure  3:  Map  of  bottom  friction  coefficients  Cf  based  on  sediment  types  (Dutkiewicz  et  al.,  2015)  with 
assumed  grain  size  and  sediment  density  (Table  2).  The  empirical  equations  (Rijn,  2007)  also  take  into 
account  depth  and  tidal  velocities. 


Firstly,  the  location  and  width  of  the  sponge  layer  l  must  be  specified.  We  take  l  to  be  approx¬ 
imately  equal  to  10%  the  length  of  the  M2  tidal  wave,  Am2  •  The  overall  solution  was  found  to  be 
fairly  insensitive  to  the  choice  of  sponge  layer  width,  but  for  l  <  0.1  Am2  the  solutions  may  not  match 
well  across  the  sponge-calculation  domain  interface.  To  show  the  location  and  width  of  the  sponge 
layer  region,  a  hatched  £+’  region  is  included  in  figures  throughout  this  paper. 

In  addition,  the  sponge  layer  requires  spatially  varying  absorption  coefficients  <t(®),  and  reference 
solutions  of  the  free  surface  r/c  and  velocities  uc.  Assuming  a  polynomial  type  function  for  the 
absorptive  coefficients  inside  the  sponge  layer,  they  are  derived  from  the  linear  shallow  water  solution: 


(T  —  (T7n 


o 


Vgh(a  +  l)\n(l/F) 

l(rc/l)a+1 


(9) 

(10) 


where  r  is  the  distance  from  the  edge  of  the  sponge  layer,  a  is  the  order  of  the  polynomial  function, 
F  is  the  reduction  factor  of  the  outgoing  wave  at  the  position  rc  from  the  edge  of  the  sponge.  The 
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parameters  a  =  2,  F  =  20  and  rc/Z  =  0.5  are  chosen  in  this  study  but  the  solution  is  not  typically 
sensitive  to  the  choice  of  these  factors.  The  reference  solutions  are  interpolated  from  TPX08  tidal 
constituents  in  the  same  way  we  specify  the  lateral  boundary.  Note  that  to  get  uCl  the  conservative 
transport  variable  ( uch )  is  interpolated  from  TPX08  before  dividing  this  by  the  nodal  depths  on 
our  grid  for  additional  consistency. 

3.7  Finite-Element  Solution 

ADCIRC  solves  the  governing  equations  in  a  continuous- Galerkin  framework,  where  the  generalized 
wave  continuity  equation  (GWCE)  is  utilized  to  eliminate  spurious  modes  (c.f.  Westerink  et  al., 
1992).  The  two-part  symmetrical  velocity  based  method  for  the  lateral  stress  terms  (Dresback, 
Kolar,  and  Luettich,  Jr.,  2005),  and  explicit  mass-lumping  mode  are  used  to  solve  the  GWCE  in 
this  study. 

A  time  step  At  =  2  s  can  be  used  with  our  current  grid  without  generating  Courant-Friedrichs- 
Lewy  (CFL)  induced  numerical  instabilities.  Wall-clock  times  are  approximately  11  min  day-1  of 
simulation  time  using  960  computational  cores  (»  10, 000  finite-element  nodes  per  core)  of  a  high- 
performance  computing  machine  with  Haswell  processors  and  a  Mellanox  FDR  Infiniband  network 
connection.  To  validate  the  model  with  observations,  we  simulated  for  195  days,  including  a  15  day 
spin-up  from  a  completely  zero  state.  The  final  180  days  are  used  for  the  harmonic  analysis  of  the 
tides.  The  long  six  month  time  period  is  required  to  correctly  separate  all  the  tidal  constituents  of 
interest  (e.g  Ki  and  Pi). 


4  Summary  of  Tidal  Validation  from  Best  Model  Setup 

4.1  Best  Model  Setup 

To  obtain  the  best  model  setup  we  first  found  the  global  amplification  factor  of  the  internal  tide 
dissipation  parameter  so  that  the  model  skill  versus  TPX08  was  maximized  in  the  deep  ocean 
(h  >  500  m).  A  positive  definite  and  spatially  smoothed  Nb  modified  version  of  the  Nycander 
method  with  C^yc  =  2.9  and  local  multiplier  coefficients  over  the  Luzon  Strait  (see  §5.3)  was 
decided  on.  Local  bathymetry  datasets  and  hand-edits  were  applied  to  shallow  water  regions  and 
responses  against  coastal  tide  gauges  were  checked  for  reliability  in  the  harmonic  analysis.  Finally,  a 
map  of  varying  bottom  friction  dissipation  coefficients  Cf  was  calculated  based  on  some  information 
of  the  local  sediment  types,  as  described  in  §3.4,  in  an  attempt  to  increase  the  model  skill  versus 
using  a  spatially  constant  Cf.  The  best  model  setup  is  denoted  by  lComp  +  IT  +  SV\  indicating 
the  use  of  our  comprehensive  bathymetric  data  (Table  1),  optimal  internal  tide  dissipation,  and 
spatially  varying  Cf. 

4.2  Measure  of  Model  Skill 

To  measure  the  skill  of  the  model  for  the  purpose  of  determining  and  evaluating  the  best  setup, 
we  compare  with  the  root- mean- square  (RMS)  discrepancy  D  of  the  elevation  (either  for  a  single 
tidal  constituent  or  for  the  total  free  surface)  at  a  point.  D  is  the  average  of  the  squared  differences 
between  measured  and  observed  elevations  integrated  over  a  long  period  of  time.  It  is  calculated 
in  this  study  using  the  sum  of  the  vector  differences  of  the  in-phase  (Afccos#fc)  and  quadrature 
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(Ak  sm6k)  components  of  each  constituent  (Wang  et  al.,  2012): 


D  =  0.5  £  [( 4 )2  +  (44  -  2 Ak0Akm  cos(0k  -  0* )] 


1/2 


(11) 


where  Ak  and  9k  are  the  amplitudes  and  phase  lags  of  the  kth  constituent  respectively,  and  the 
subscripts  ‘o’  and  ‘m’  refer  to  the  observed  and  modeled  values  respectively.  In  addition,  the  relative 
RMS  discrepancy  RD  =  D/V,  where  V  is  the  absolute  average  value  of  the  variability  in  the  free 
surface  elevation,  and  is  calculated  by  (Wang  et  al.,  2012): 


V  = 


o.5£(4)2 


(12) 


For  an  overview  of  the  spatial  distribution,  we  include  scatter  plots  of  D  and  RD  at  the  tide 
gauges  (and  contour  plots  versus  TPX08)  in  order  to  highlight  regions  of  notably  small  or  large 
discrepancies.  However,  to  obtain  a  single  global  metric  of  performance  the  mean  of  the  discrepancy 
D,  denoted  D ,  or  the  mean  of  RD ,  denoted  RD ,  is  used.  Note  that  when  calculating  D  over  a 
region  to  compare  against  TPX08  ( tpx )  this  is  computed  as: 


D 


JfDdA 

tpx~-ffdA 


(13) 


where  ff  dA  indicates  an  area  integral  that  is  performed  over  the  elements  of  the  grid.  When 
comparing  against  tide  gauges  ( Dtg ,  RDtg),  the  arithmetic  average  is  used.  In  comparison  to  D, 
the  RMSE  metric  commonly  used  (e.g  Stammer  et  al.,  2014;  Buijsman  et  al.,  2015)  is: 


RMSE=/^lr  (14) 

i.e.  it  is  the  square-root  of  the  mean  of  D2  and  is  always  larger  than  D.  The  RMSE  may  experience 
abrupt  changes  with  depth  and  tends  to  overestimate  the  overall  discrepancy  (Wang  et  al.,  2012). 
In  contrast,  D  has  been  shown  to  decrease  monotonically  with  depth  (Wang  et  al.,  2012),  thus  we 
chose  to  predominantly  use  D.  However,  we  also  quote  values  of  RMSE  for  comparison  with  those 
reported  in  other  studies. 

4.3  Tidal  Gauge  Database 

A  database  of  harmonic  constituents  consisting  of  39  deep-water  stations,  62  shallow  water/shelf 
stations,  and  659  unique  coastal  tide  gauge  locations  was  assembled  from  multiple  sources  for  the 
computational  domain  (detailed  in  Table  3).  Within  the  coastal  tide  gauge  set  there  are  a  number 
of  data-points  duplicated  between  sources  so  we  set  up  a  hierarchy  between  the  different  sources  to 
decide  what  value  to  use  in  our  model  evaluation  based  on  perceived  reliability  (Table  3  is  listed  in 
hierarchical  order,  and  the  number  of  stations  listed  for  each  source  is  the  eventual  number  after 
removal  of  duplicates). 

The  KHOA  (Korean  Hydrographic  and  Oceanographic  Agency),  GESLA-2  (Global  Extreme  Sea 
Level  Analysis),  and  UHSLC  FD  (University  of  Hawaii  Sea  Level  Center  fast- delivery)  databases 
contain  long-term  hourly  time  series  of  elevations.  We  have  used  the  Utide  MATLAB  function 


11 


utsolv  (Codiga,  2011),  which  uses  the  iteratively- weighted  least-square  harmonic  analysis  tech¬ 
nique,  to  obtain  up  to  68  tidal  constituents.  Within  some  of  the  sources,  not  all  of  the  main  eight 
constituents  are  necessarily  posted  such  as  in  the  JMA,  NBoB,  SCS,  Yellow  Sea,  most  of  IHO,  and 
TOPEX/POSEIDON  Crossover  sets  (M2,  S2,  Ki,  Oi  only).  The  Java  Sea/SCS  set  excludes  K2  and 
Pi.  The  ST727  excludes  various  constituents  (often  just  Qi)  at  approximately  half  of  the  locations. 
Occasionally,  locations  in  the  GESLA.2,  and  Truth-Shallow  sets  miss  one  or  two  of  the  eight  main 
constituents  due  to  short  analyses.  Additionally,  posted  geographical  locations  are  often  to  the 
nearest  1'  or  0.1',  so  we  have  checked  station  positions  to  ensure  that  they  are  likely  to  be  in  the 
correct  location,  in  particular  not  on  land.  The  FUGAWI  navigational  charts  often  indicate  known 
locations  of  tide  gauges  and  were  used  to  correct  locations  where  applicable.  Some  stations  were 
excluded  because  they  were  located  just  outside  of  our  computational  grid  (in  back-bays,  rivers, 
etc.).  The  locations  and  constituent  values  used  to  conduct  the  statistical  comparisons  in  this  study 
(sans  proprietary  data)  can  be  found  at  Pringle  (2017). 

A  detailed  discussion  on  the  harmonic  analyses  and  uncertainties  of  the  tidal  constituents  (ex¬ 
tended  globally)  will  be  reported  in  a  subsequent  paper.  For  example,  the  constituents  from  the 
IHO  and  ST727  gauges  are  derived  between  1848-1970  and  may  have  larger  error  margins.  In  addi¬ 
tion,  the  rate  at  which  the  coasts  are  changing  globally  is  accelerating,  and  elimination  of  wetlands, 
land  reclamation,  and  deep  dredging  of  channels  and  ports  may  impact  tidal  constituents  which  are 
typically  located  within  harbor  complexes  in  support  of  navigation. 

4.4  Overall  Global  Assessment 

We  start  by  showing  the  global  responses  of  the  M2  and  Ki  (Fig  4)  tidal  waves,  and  their  RMS 
discrepancies  against  TPX08  ( Dtpx )  for  the  Comp  +  IT  +  SV  model  setup.  The  general  response  is 
well  described  for  both  constituents  by  our  model  including  positions  of  most  amphidromes.  Obvious 
exceptions  are  the  M2  amphidromes  off  the  south-west  tip  of  Australia  and  near  Mawson  Station, 
Antarctica.  The  positions  of  these  amphidromes  and  the  solution  in  the  Southern  Ocean  were  found 


Table  3:  Tide  gauge  data  sources,  number  and  availability.  Listed  in  hierarchical  order  for  the  coastal 
gauges 


Name 

Source 

Number 

Type 

Availability 

Truth-Pelagic 

Shum  et  al.  (1997) 

31 

deep-water  const. 

free  at  website1 

Truth-Shallow 

Stammer  et  al.  (2014) 

52 

shallow-water  const. 

free  at  website1 

TOPEX/POSEIDON  Crossovers 

Robertson  and  Ffield  (2008) 

8/5 

deep/shallow-water  const. 

listed  in  paper 

Java  Sea/SCS 

Wei  et  al.  (2016) 

5 

shallow-water  const. 

listed  in  paper 

NOAA 

NOAA  Tides  and  Currents  (2017) 

4 

coastal  const. 

free  at  website2 

JMA 

Japanese  Meteorological  Agency  (2017) 

181 

coastal  const. 

free  at  website3 

AusTides 

Australian  National  Tide  Tables  (2013) 

63 

coastal  const. 

proprietary 

KHOA 

Korean  Hydrographic  and  Oceanographic  Agency  (2017)  35 

coastal  elev. 

free  at  website4 

GESLA-2 

Woodworth  et  al.  (2017) 

107 

coastal  elev. 

free  at  website5 

UHSLC  FD 

Caldwell,  Mcrriiicld.  and  Thompson  (2015) 

19 

coastal  elev. 

free  at  website6 

NBoB 

Krien  et  al.,  2016 

2 

coastal  const. 

listed  in  paper 

SCS 

Fang  et  al.  (1999) 

29 

coastal  const. 

listed  in  paper 

Yellow  Sea 

Fang  et  al.  (2004) 

6 

coastal  const. 

listed  in  paper 

ST727 

British  Hydrographic  Institute  (c. 1848-1970) 

130 

coastal  const. 

free  at  website1 

IHO 

International  Hydrographic  Office  (1990) 

83 

coastal  const. 

proprietary 

const.:  indicates  data  is  available  directly  as  harmonic  constituents 
elev.:  indicates  data  is  available  as  time  series  of  elevations 

1 :  ftp : //ftp . legos . obs-mip . f r/pub/FES2012-pro ject/data/gauges/2013- 12- 16/ 
2 :  https : / /t idesandcurrents . noaa . gov/ gmap3/ 

3:  http : / /www . data . jma . go . jp/kaiyou/db/tide/suisan/station2017 . php 
4:  http : //www . khoa . go . kr/koof s/kor/observat ion/obs_real . do 
5:  http://www.gesla.org/ 

6:  ftp://ftp.soest.hawaii.edu/uhslc/fast 


12 


Figure  4:  Amplitude  (m)  and  phase  responses  of  the  (i)  M2  and  (ii)  Ki  tidal  waves;  (a)  Comp  +  IT  +  SV 
model  setup,  (b)  TPX08  model  (obtained  from  http://volkov.oce.orst.edu/tides/tpxo8_atlas.html), 
(c)  RMS  discrepancies  (m)  between  Comp  -h  IT+  SV  model  setup  and  TPX08,  Dtpx.  ‘-f5  hatched  regions 
indicate  absorptive  sponge  zone. 


13 


Figure  5:  Spatial  distribution  of  discrepancies  of  the  (i)  M2  and  (ii)  Ki  tidal  waves  versus  tide  gauges 
for  the  Comp  +  IT  +  SV  model  setup;  (a)  RMS  discrepancy  Dtg ,  (b)  relative  RMS  discrepancy  RDtg . 
Triangles:  deep  water  gauges,  Squares:  continental  shelf  water  gauges,  Circles:  coastal  tide  gauges. 


to  be  very  sensitive  to  the  boundary  conditions  applied  in  this  study  and  may  be  impacted  by  the 
TPX08  derived  velocities  in  the  absorption-generation  sponge  layer  (see  §5.1). 

The  spatial  distribution  of  the  RMS  discrepancies  ( Dtg )  and  relative  discrepancies  ( RDtg )  for  the 
Comp  +  IT  +  SV  model  setup  against  tide  gauges  are  also  illustrated  (Fig  5).  Overall,  tide  gauges 
with  similar  discrepancies  are  generally  clustered  together,  and  there  is  a  relatively  strong  spatial 
correlation  between  discrepancies  against  TPX08  and  those  at  tide  gauges.  Exceptions  to  this 
include  much  of  the  inner  coast  of  the  Yellow  Sea  and  the  Seto  Inland  Sea  where  the  TPX08  model 
may  not  be  reliable.  The  Comp  +  IT  +  SV  model  setup  performs  particularly  well  throughout 
the  western  Pacific  Ocean  including  along  the  Japanese  archipelago  and  northeastern  Australia  for 
both  constituents.  Notable  wide-spread  RMS  discrepancies  in  the  M2  tidal  wave  appear  in  the 
Mozambique  Channel,  north  and  west  Arabian  Sea,  Red  Sea,  Sea  of  Okhotsk,  Andaman  Sea,  Yellow 
Sea,  northern  Australian  shelf  and  the  Celebes  Sea.  Ki  RMS  discrepancies  are  notable  in  the  Sea 
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of  Okhotsk,  Arabian  Sea,  South  China  Sea  and  Java  Seas,  and  the  Arafura  Sea.  Predominantly 
large  tidal  ranges  account  for  the  discrepancies  shown.  For  example,  M2  RDtg  values  are  relatively 
small  in  the  Yellow  Sea  even  though  Dtpx  values  appear  large  in  the  Yellow  Sea  for  the  Comp  +  IT 
+  SV  model  setup.  In  fact,  the  response  is  improved  rather  substantially  from  the  Comp  +  IT  -f 
SC  model  setup  here  (§5.4).  RDtg  is  also  less  significant  than  Dtg  in  the  Mozambique  Channel  and 
northern  Australian  shelf.  These  two  regions  are  heavily  influenced  by  large-scale  effects  related  to 
lateral  boundary  conditions  (§5.1)  and  internal  tide  dissipation  (§5.3). 

On  the  other  hand,  both  Dtg  and  RDtg  are  large  in  the  Sea  of  Okhotsk  for  both  constituents. 
The  importance  of  bathymetry  in  the  region  (which  is  not  well  known)  has  been  highlighted  by 
Zaron  (2017).  The  Celebes  Sea  is  also  a  problem  area  for  M2  that  is  most  likely  a  result  of  incorrect 
flux  exchanges  through  the  island  chains  due  to  inadequate  bathymetry  and  a  poor  representation 
of  internal  tide  dissipation  particularly  in  shallow  waters.  It  should  be  noted  that  the  Celebes  Sea 
and  surrounding  Indonesian  seas  was  a  focus  of  the  original  TPXO  study  (Egbert  and  Erofeeva, 
2002)  due  to  its  poor  forward  model  responses,  and  the  region  has  been  found  to  cause  problems 
for  three-dimensional  ocean  circulation  models  (Robertson  and  Ffield,  2008;  Ngodock  et  al.,  2016). 
The  South  China  and  Java  Sea  region  extending  down  to  the  Torres  Strait  has  a  relatively  large 
diurnal  tidal  range  and  Ki  Dtg  and  RDtg  values  are  not  small  compared  to  most  of  the  domain. 
The  physics  of  the  Ki  tidal  wave  here  can  be  thought  of  as  a  standing  wave  where  the  response  is 
likely  to  depend  highly  on  the  overall  bathymetry  and  shoreline  of  the  region.  The  region  is  also 


Figure  6:  Amplitudes,  A  (left)  and  phase  lags,  0  (right)  of  up  to  all  eight  major  tidal  constituents  for  the 
Comp  +  IT  +  SV  model  setup  versus  observed  values  at  tide  gauges.  Triangles:  deep  water  gauges,  Squares: 
continental  shelf  water  gauges,  Circles:  coastal  tide  gauges.  Colors  of  markers  for  the  amplitude  refer  to 
the  absolute  error  (m)  between  model  and  observed.  Colors  of  markers  for  the  phase  refer  to  the  absolute 
errors  normalized  by  180°  between  model  and  observed.  The  same  color  scale  as  Fig  5(a)  is  used  for  both. 
Statistics  shown  on  the  figure  are  as  follows:  R2  is  the  coefficient  of  determination,  ostd  is  the  standard 
deviation  of  the  error,  E  is  the  mean  error,  \E\  is  the  mean  absolute  error,  and  En  is  the  normalized  mean 
absolute  error. 
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heavily  influenced  by  the  energy  flux  permitted  through  the  Luzon  Strait  which  is  largely  controlled 
by  internal  tide  dissipation  (§5.3).  Note  that  in  some  areas  such  as  between  South  China  Sea  and 
Java  Sea,  which  has  a  small  M2  tidal  range  because  it  is  close  to  an  amphidrome,  RDtg  becomes 
very  large  however  Dtg  is  relatively  small. 

A  summary  of  the  global  tide  gauge  errors  shown  in  terms  of  amplitudes  ( R 2  =  0.93,  <jstd  —  0.09 
m,  \E\  =  0.04  m)  and  phases  (R2  =  0.97,  <jstd  =  18.3°,  \E\  =  10.3°)  of  up  to  all  eight  major  tidal 
constituents  from  the  Comp  +  IT  +  SV  model  setup  against  the  observed  values  is  presented  (Fig  6, 
see  caption  for  definitions  of  error  metrics).  There  are  a  total  of  6080  data  points  on  each  plot.  Just 
2.4%  of  them  represent  absolute  amplitude  errors  >  0.2  m,  and  2.9%  represent  absolute  phase  errors 
>  36°  (colored  orange  to  purple  in  Fig  6).  Outliers  in  the  amplitudes  of  constituents  tend  to  be 
underestimates  rather  than  overestimates  which  may  indicate  deltaic  regions,  estuaries,  back-bays 
and  rivers  where  the  bathymetry  is  inadequate  and  overly  dissipative,  e.g.  the  Ganges  Delta  where 
large  discrepancies  are  present  (Fig  5).  Aside  from  these  regions,  there  is  a  consistent  spread  of 
errors  for  both  the  amplitudes  ( E  =  -0.02  m)  and  phases  (E  =  1.37°)  indicating  a  largely  unbiased 
system. 

4.5  Assessment  in  Deep,  Shelf,  and  Coastal  Waters 

This  section  summarizes  the  results  by  dividing  them  into  regions  based  on  depth  and  proximity  to 
the  coastline.  The  three  regions  are  as  follows:  deep  water  ( h  >  500  m),  continental  shelf  and  slope 
waters  (25  <  h  <  500  m),  and  coastal  waters  (includes  continental  and  island  coastlines).  Table  4 
compares  discrepancies  between  the  various  IndWPac  model  setups  (different  bathymetry  datasets, 
with  and  without  internal  tide  dissipation  and  spatially  varying  bottom  friction  coefficients),  at 
the  tide  gauges  and  the  TPX08  atlas  for  these  three  regions.  Note  that  when  the  interpolation 
of  TPX08  to  the  coastal  tide  gauges  was  performed  using  their  native  data  extraction  program 
OTPS2,  a  total  of  93  locations  returned  a  null  value.  Thus,  for  a  fair  comparison  we  present  our 
model  results  against  this  reduced  set  of  stations.  The  statistics  of  the  IndWPac  model  were  not 
noticeably  different  for  the  full  coastal  gauge  set. 

4.5.1  Discrepancy  in  Deep  Water 

The  total  free  surface  mean  discrepancies  for  the  Comp  +  IT  +  SV  model  setup  (Dtg  =  4.7  cm, 
RDtg  =  13%)  are  2.4  times  those  of  the  TPX08  atlas  (Dtg  =  2.0  cm,  RDtg  =  5.5%)  at  the  deep 
water  tide  gauges.  The  IndWPac  model  discrepancies  in  deep  water  are  predominantly  affected  by 
the  internal  tide  dissipation  which  reduces  the  total  free  surface  RMS  discrepancy  by  47%  ( Dtpx )  and 
55%  ( Dtg ).  Different  bathymetry  datasets  and  bottom  friction  coefficients  have  little  effect.  Hot¬ 
spots  of  discrepancy  against  deep  water  tide  gauges  for  the  IndWPac  model  occur  in  the  Celebes  Sea 
and  Banda  Sea  (see  Fig  5)  against  TOPEX/POSEIDON  satellite  crossover  observations  (Robertson 
and  Ffield,  2008),  particularly  for  the  M2  tidal  wave.  Without  the  crossover  points  (which  are 
technically  not  tide  gauges)  the  M2  Dtg  for  the  Comp  -h  IT  -f  SV  model  setup  is  closer  to  2  cm 
instead  of  3.6  cm. 

A  non-assimilative  hydrodynamic  model  (Buijsman  et  al.,  2015)  achieves  an  RMSE  of  approx¬ 
imately  4  cm  in  the  Indian  and  Pacific  Oceans,  and  most  hydrodynamics  models  give  an  RMSE 
greater  than  5  cm  (Stammer  et  al.,  2014)  globally  for  the  M2  tidal  wave  in  waters  deeper  than  1000 
m  against  TPX08.  In  comparison,  the  Comp  -h  IT  +  SV  model  setup  gives  RMSE  =  3.6  cm. 
While  the  discrepancy  for  our  model  is  somewhat  smaller  than  other  hydrodynamic  models,  it  is  not 
a  significant  change  despite  generally  higher  resolution  of  the  grid  and  nearshore  bathymetric  data. 
As  shown,  only  internal  tide  dissipation  resulted  in  a  notable  reduction  to  the  deep  water  discrepan- 
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Table  4:  The  mean  RMS  ( Dtg )  and  relative  RMS  ( RDtg )  discrepancies  of  the  M2,  Ki,  and  the  total  free 
surface  (up  to  all  eight  major  constituents  combined)  at  tide  gauges  for  various  IndWPac  model  setups  plus 
the  TPX08  atlas  (http://volkov.oce.orst.edu/tides/tpxo8_atlas.html),  separated  into  three  different 
regions  (deep,  continental  shelf  and  slope,  and  coastal).  The  mean  RMS  discrepancy  against  TPX08  ( Dtpx ) 
is  shown  in  deep,  and  continental  shelf  and  slope  waters.  Stations  numbers,  units,  and  standard  deviations 
are  in  parentheses 


Model 

Error 

Tidal 

Comp  -h 

GEBCO  + 

SRTM  -h 

Comp  + 

Comp  + 

TPYOS 

Metric 

Wave 

NoIT  -h  SC 

IT  -h  SC 

IT  +  SC 

IT  -f  SC 

IT  +  SV 

J.  1  AvJO 

m2 

5.69  (4.31) 

2.97  (2.21) 

2.92  (2.26) 

2.90  (2.20) 

2.89  (2.18) 

- 

Dtpx 

Ki 

1.39  (1.56) 

1.08  (3.14) 

1.01  (3.59) 

0.95  (1.19) 

0.95  (1.17) 

- 

(cm) 

All 

6.91  (4.76) 

3.85  (5.51) 

3.82  (6.23) 

3.67  (2.59) 

3.67  (2.56) 

- 

m2 

8.90  (9.97) 

3.90  (4.40) 

3.82  (4.02) 

3.66  (4.45) 

3.55  (4.34) 

0.86  (0.87) 

Deep 

Dtg 

Ki 

2.16  (1.66) 

1.03  (0.82) 

0.94  (0.71) 

0.93  (0.71) 

0.92  (0.68) 

0.50  (0.34) 

(39) 

(cm) 

All 

10.8  (11.0) 

5.09  (5.03) 

5.02  (4.60) 

4.82  (4.98) 

4.71  (4.89) 

2.02  (2.82) 

m2 

32.7  (25.7) 

15.7  (11.4) 

15.7  (11.3) 

14.2  (11.9) 

13.9  (11.7) 

3.82  (3.31) 

RDtg 

Ki 

16.4  (9.08) 

8.40  (6.31) 

7.72  (5.72) 

7.67  (5.84) 

7.69  (5.83) 

4.29  (3.34) 

(%) 

All 

28.4  (19.8) 

14.3  (10.2) 

14.3  (9.35) 

13.4  (10.3) 

13.1  (10.1) 

5.46  (5.67) 

m2 

11.2  (11.0) 

8.05  (9.02) 

7.58  (9.45) 

6.76  (7.67) 

6.48  (7.76) 

- 

Dtpx 

Ki 

5.83  (13.7) 

6.69  (14.4) 

6.57  (17.2) 

4.52  (7.08) 

4.75  (7.71) 

- 

(cm) 

All 

16.6  (22.1) 

15.0  (23.6) 

14.8  (28.7) 

10.9  (10.8) 

11.0  (11.4) 

- 

m2 

18.4  (12.1) 

12.3  (11.4) 

14.9  (19.3) 

9.70  (9.63) 

9.35  (9.87) 

2.91  (3.28) 

Shelf 

Dtg 

Ki 

5.71  (4.54) 

4.41  (4.05) 

4.84  (5.04) 

3.71  (4.13) 

4.47  (4.78) 

1.60  (1.56) 

(62) 

(cm) 

All 

22.8  (12.3) 

16.0  (11.7) 

19.4  (20.8) 

13.0  (11.0) 

13.4  (11.4) 

5.41  (3.76) 

m2 

76.9  (89.1) 

53.8  (73.3) 

52.6  (65.4) 

42.1  (58.5) 

40.9  (55.8) 

12.8  (20.9) 

RDtg 

Ki 

24.4  (14.5) 

21.0  (18.9) 

21.5  (20.0) 

19.0  (20.0) 

19.9  (20.0) 

7.89  (9.37) 

(%) 

All 

40.2  (18.1) 

26.8  (15.0) 

29.9  (19.8) 

22.1  (14.8) 

22.6  (16.1) 

9.24  (5.37) 

m2 

16.0  (17.1) 

24.6  (33.4) 

17.9  (24.1) 

12.1  (15.6) 

10.5  (14.4) 

13.5  (39.3) 

Dtg 

Ki 

4.58  (6.15) 

6.95  (9.18) 

5.45  (7.80) 

4.09  (6.34) 

3.90  (6.32) 

3.12  (5.58) 

Coast 

(cm) 

All 

20.9  (19.6) 

29.9  (37.1) 

22.5  (27.4) 

15.8  (18.4) 

14.4  (17.2) 

17.0  (47.9) 

(659) 

m2 

45.9  (60.4) 

45.7  (44.0) 

36.3  (40.1) 

28.5  (36.9) 

27.2  (38.4) 

24.7  (38.6) 

RDtg 

Ki 

27.0  (38.2) 

37.9  (33.8) 

29.0  (27.4) 

22.4  (19.9) 

21.6  (20.5) 

18.1  (21.4) 

(%) 

All 

36.4  (22.5) 

41.6  (31.9) 

32.8  (26.0) 

25.3  (18.2) 

24.4  (19.2) 

22.8  (30.8) 

Model  Setups 

Bathymetry:  ‘ GEBCO ’  uses  GEBCO_2014  bathymetric  data,  ‘ SRTM ’  uses  SRTM15JPLUS  bathymetric  data, 
lCompy  uses  our  comprehensive  bathymetric  data  (Table  1) 

Internal  Tide  Dissipation:  'NoITy  does  not  include  internal  tide  dissipation,  lITy  uses  the  optimal  internal  tide 
parameters 

Bottom  Friction:  lSCy  uses  a  spatially  constant  Cf  —  2.5xl0-3,  lSV*  uses  the  spatially  varying  Cf  map  (Fig.3) 
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cies.  Better  nearshore  bathymetry  and  grid  resolutions  do  not  allow  for  significant  improvements  in 
the  internal  tide  dissipation  matrix  compared  with  coarser  grid  models  because  its  calculation  relies 
mostly  on  the  deep  water  satellite  altimetry  data  in  global  bathymetric  datasets  that  is  still  limited 
to  >  10  km  resolution  accuracy  (Goff  and  Arbic,  2010).  Furthermore,  the  topographic  roughness 
can  be  calculated  on  the  relatively  fine  ~1  km  bathymetric  grid  before  being  interpolated  onto  the 
coarser  computational  grid  for  parameterization  in  barotropic  models,  reducing  the  requirement  for 
a  fine  grid  in  the  ocean. 

4.5.2  Discrepancy  in  Continental  Shelf  and  Slope  Waters 

The  total  free  surface  Dtg  at  tide  gauges  on  the  continental  shelf  are  2.6  to  2.8  times  larger  than 
those  in  deep  water  for  both  the  Comp  +  IT  +  SV  model  setup  and  the  TPX08  atlas.  Similar  to 
deep  water  regions,  the  Comp  +  IT  -h  SV  model  total  free  surface  discrepancies  ( Dtg  =  13  cm,  RDtg 
—  23%)  are  2.5  times  those  of  the  TPX08  atlas  ( Dtg  =  5.4  cm,  RDtg  =  9.2%)  at  the  tide  gauges 
on  the  continental  shelf.  The  most  significant  factors  in  reducing  the  total  free  surface  discrepancies 
are  internal  tide  (34%  reduction  in  Dtpx)  and  the  nearshore  bathymetric  datasets  (26%  reduction  in 
Dtpx).  The  bottom  friction  coefficient  has  a  small  impact  overall,  although  the  discrepancy  for  the 
Ki  constituent  increased  when  using  the  spatially  varying  Cj  map.  The  GEBCO-2014  bathymetry 
model  gives  lower  Dtg  than  the  SRTM15.PLUS  model,  however  Dtpx  are  quite  similar  between  the 
two  bathymetric  datasets.  Note  that  Dtpx  tends  to  give  a  smoother  indicator  of  the  change  between 
model  setups  because  it  is  integrated  over  the  whole  domain  (where  25  <  h  <  500  m).  Furthermore, 
Dtpx  is  2.9  cm  or  ~30%  smaller  than  Dtg  for  M2.  This  could  be  because  the  shelf  gauges  tend  to 
be  in  regions  with  large  tidal  ranges  such  as  the  northern  regions  of  the  Australian  shelf  and  the 
Yellow  Sea. 

In  comparison  to  other  forward  hydrodynamic  models  the  performance  is  favorable.  Stammer  et 
al.  (2014)  reports  that  the  global  M2  RMSE  in  shelf  waters  is  24-49  cm  against  tide  gauges  and  19-28 
cm  versus  TPX08.  Comparatively,  the  M2  RMSE  is  13.1  cm  against  tide  gauges  and  10.1  cm  versus 
TPX08  for  the  Comp  -h  IT  +  SV  model  setup.  Although  it  should  be  kept  in  mind  that  the  Stammer 
et  al.  (2014)  errors  are  global  and  hence  they  cannot  be  treated  as  a  direct  comparison,  according 
to  the  TPX08  atlas  the  total  energy  density  ( TED ,  defined  in  §5.1)  of  the  M2  tidal  wave  is  slightly 
higher  in  the  IndWPac  domain  (820  Jm“2)  compared  to  the  entire  globe  (695  Jm-2).  So  while  data- 
assimilation  clearly  provides  greater  accuracy  on  the  continental  shelf  than  forward  hydrodynamic 
models  can  currently  manage,  computational  grids  that  become  fine  enough  to  adequately  support 
high-resolution  bathymetric  data  nearshore  and  on  the  shelf  (as  in  the  IndWPac  model)  do  appear 
to  result  in  noticeably  smaller  errors  than  coarser  forward  hydrodynamic  models. 

4.5.3  Discrepancy  at  the  Coast 

The  RMS  discrepancies  at  the  coastal  tide  gauges  are  only  marginally  larger  than  those  on  the  shelf 
for  the  Comp  +  IT  +  SV  model  setup.  However,  the  discrepancies  increase  significantly  from  the 
shelf  to  the  coast  when  only  the  global  bathymetry  datasets  are  used,  in  particular  GEBCCL2014. 
The  nearshore  bathymetry  dataset  plays  a  large  role  in  reducing  the  discrepancy  (30%  reduction  in 
total  free  surface  Dtg).  The  spatially  varying  Cj  map  has  a  smaller  but  noticeable  global  effect  (9% 
reduction  in  total  free  surface  Dig).  Local  effects  of  Cj  are  detailed  in  §5.4.  At  approximately  77% 
of  the  tide  gauge  locations  the  total  free  surface  Dtg  of  the  Comp  +  IT  +  SV  model  setup  is  less 
than  20  cm  (Fig  7).  This  is  the  same  as  the  target  metric  used  for  a  high-resolution  western  North 
Atlantic  model  (much  smaller  in  scale  to  IndWPac)  where  this  was  satisfied  at  324  of  398  (81%) 
locations  (Technology  Riverside  Inc.  and  AECOM,  2015). 
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Figure  7:  Cumulative  distribution  functions  of  the  total  free  surface  (up  to  all  eight  tidal  constituents) 
RMS  discrepancies  Dtg  (left),  and  relative  RMS  discrepancies  RDtg  (right),  versus  coastal  tide  gauges  for 
different  IndWPac  model  setups  (see  explanation  in  footnotes  of  Table  4)  and  the  TPX08  atlas. 


Even  though  the  total  free  surface  Dtg  for  the  Comp  +  IT  +  SV  model  setup  is  2.5  times  that 
of  the  TPX08  atlas  on  the  shelf,  Dtg  is  2.6  cm  (15%)  smaller  than  the  TPX08  atlas  for  the  Comp 
+  IT  +  SV  model  setup  at  the  coast.  However,  as  illustrated  by  the  cumulative  distribution  (cdf) 
curves  (Fig  7),  a  higher  percentage  of  locations  will  be  within  a  given  target  discrepancy  up  to  Dtg 
—  24  cm  ( RDtg  =  35%)  for  the  TPX08  model.  On  the  other  hand,  the  TPX08  model  cdf  curves 
have  long  tails  indicating  a  number  of  high-magnitude  outliers,  whereas  this  is  not  the  case  for  the 
IndWPac  model  with  the  nearshore  bathymetry  included.  A  related  trend  is  also  evident:  while  Dtg 
for  the  Comp  +  IT  -h  SV  model  setup  model  at  the  coast  is  only  marginally  increased  from  those 
on  the  shelf  (and  even  decreased  for  Ki),  Dtg  for  the  TPX08  atlas  has  increased  almost  two-fold  for 
Ki  and  five-fold  for  M2.  Furthermore,  the  standard  deviations  of  Dtg  for  TPX08  are  much  larger 
than  Dtg.  In  contrast,  the  standard  deviations  of  Dtg  for  the  IndWPac  model  setups  using  our 
comprehensive  bathymetric  data  are  similar  in  magnitude  to  Dtg . 

The  total  free  surface  RDtg  are  slightly  smaller  for  the  TPX08  atlas  (23%)  than  the  Comp  + 
IT  +  SV  model  setup  (24%).  Furthermore,  the  IndWPac  model  cdf  curves  for  RDtg  start  off  flat 
and  begin  from  a  higher  minimum  value  than  those  of  the  TPX08  atlas,  whereas  this  is  not  the 
case  for  the  Dtg  curves.  This  could  be  explained  by  the  difficulty  in  obtaining  accurate  relative 
discrepancies  in  regions  of  small  tides  especially  those  near  amphidromes,  e.g.  in  southern  Australia 
RDtg  becomes  large  due  to  misplacement  of  the  amphidromc  off  the  southwest  tip  of  Australia  even 
though  Dtg  remain  somewhat  small  (Fig  5).  Free  surface  elevations  associated  with  amphidromes 
that  are  situated  offshore  like  this  can  be  better  placed  through  data-assimilation  by  e.g.  TPX08. 

Overall,  the  large  numbers  and  magnitudes  of  outliers  and  standard  deviations  of  RMS  discrep¬ 
ancies  likely  indicate  that  the  TPX08  atlas  is  unable  to  capture  small-scale  changes  in  amplitude 
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and/or  phase  that  can  occur  in  bays  and  harbors  or  in-behind  small  islands  and  peninsulas  due  to 
coarse  resolution  and  bathymetry.  But  since  TPX08  cdf  curves  rise  quickly  initially,  there  is  an 
indication  that  in  coastal  regions  where  the  solution  is  not  significantly  different  from  that  offshore 
(and  where  the  gauges  have  been  included  in  the  assimilation  process)  TPX08  is  accurate.  How¬ 
ever,  although  TPX08  use  1/30°  resolution  local  insets  to  improve  the  coastal  solutions,  even  finer 
computational  grids  on  the  coast  that  support  high-resolution  local  bathymetry  may  be  required  to 
make  further  inroads  into  coastal  performance  elsewhere,  particularly  where  nonlinearities  play  an 
important  role  (c.f.  Lefevre,  Provost,  and  Lyard,  2000;  Lyard  et  al.,  2006). 

In  contrast,  the  statistics  indicate  that  the  high-resolution  computational  grids  and  bathymetric 
data  included  in  our  model  leads  to  the  ability  to  better  capture  the  faster  changing  character¬ 
istics  of  tides  (particularly  semi-diurnal  ones)  hence  there  are  fewer  large  magnitude  outliers  and 
a  smaller  mean  discrepancy  compared  with  the  TPX08  atlas.  However,  at  a  high  percentage  of 
locations,  TPX08  gives  smaller  discrepancies  because  the  IndWPac  model  suffers  from  a  larger 
residual  discrepancy  coming  from  the  offshore  dynamics  due  to  the  absence  of  data-assimilation. 
If  data-assimilation  is  not  involved,  then  it  would  appear  that  significant  improvements  in  offshore 
bathymetric  data  and  internal  tide  dissipation  dynamics  need  to  occur  to  elevate  the  median  perfor¬ 
mance  at  the  coast  in  large-scale  models.  The  following  sections  describe  the  sensitivities  to  these 
and  other  factors. 


5  Sensitivities  to  Lateral  Boundary  Conditions,  Bathymetry, 
and  Dissipative  Controls 

5.1  Lateral  Boundary  Conditions 

In  this  study  it  was  found  that  the  most  dramatic  effect  on  the  solution  occurred  when  modifying  the 
position  and/or  lateral  boundary  condition  type  (Fig.  8).  It  turns  out  that  the  boundaries  used  in 
the  final  IndWPac  model  (which  we  call  the  ‘two  open  boundaries  domain’  in  this  section)  are  fairly 
well  placed.  In  the  initial  stages  of  the  IndWPac  model,  the  domain  was  set  up  so  that  the  Western 
Pacific  boundary  was  split  into  two  separate  boundaries  where  one  of  the  boundaries  was  defined 
spanning  from  the  Great  Australian  Bight  down  to  Antarctica  parallel  with  latitude  (designated  as 
the  ‘three  open  boundaries  domain’).  To  evaluate  the  effect  of  the  lateral  boundary  we  quantify  the 
total  energy  density,  TED  of  the  kth  tidal  wave: 


TEDk 


Po. If  (h\Uk\2  +  g{Ak?)  dA 
2  JS  dA 


(15) 


where  po  is  the  reference  density  of  sea  water.  A  total  of  four  simulations  are  conducted:  the  two 
open  boundaries  domain  with  and  without  an  absorption-generation  sponge  layer,  and  the  three 
open  boundaries  domain  with  and  without  a  sponge  layer.  Note  that  none  of  these  simulations  use 
internal  tide  dissipation  and  a  spatially  constant  Cf  =  2.5xl0“3  is  employed. 


5.1.1  Boundary  Placement  Effects  without  Sponge  Layer 

For  simulations  without  the  sponge  layer,  TEDM2  =  899  Jm~2  for  the  two  open  boundaries  domain 
and  3210  Jm“2  when  using  the  three  open  boundaries  domain,  representing  a  247%  increase.  For 
comparison,  in  the  TPX08  model  TEDU2  =  820  Jm-2  within  the  two  open  boundaries  domain. 
The  79  Jm-2  excess  in  total  energy  density  between  the  two  open  boundaries  domain  solution 
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(a.ii)  M2:  Two  open  boundaries 
+  Sponge 


(a.i)  M2:  Two  open  boundaries 


(bJ)  M2  Three  open  boundaries 


(b.ii)  Mj:  Three  open  boundaries 
-1-  Sponge 


Figure  8:  Responses  of  the  M2  tidal  wave  with  no  internal  tide  dissipation,  and  spatially  constant  Cj 
=  2.5x10  3 ;  (i)  no  absorption-generation  sponge  layer,  (ii)  with  absorption-generation  sponge  layer  (‘-b’ 
hatched  regions  indicate  sponge  layer),  (a)  two  open  boundaries  domain,  (b)  three  open  boundaries  domain. 


21 


and  TPX08  can  be  explained,  in  part,  by  the  absence  of  internal  tide  dissipation.  The  effect 
on  the  solution  due  to  boundary  placement  was  not  nearly  as  prominent  for  the  diurnal  Ki  tidal 
wave.  TEDKl  =  288  Jm-2  for  the  two  open  boundaries  domain  and  298  Jm-2  for  the  three  open 
boundaries  domain,  representing  a  small  3.5%  increase.  Most  other  constituents  also  showed  single 
digit  percent  increases  except  for  the  other  two  lunar  semi-diurnal  tides,  K2  and  N2.  In  particular, 
K2  whose  response  is  most  similar  looking  to  M2,  saw  a  231%  increase.  N2  was  increased  by  31.4%. 

5.1.2  Effects  of  Sponge  Layer 

It  was  found  that  applying  the  absorption-generation  sponge  layer  even  for  the  three  open  boundaries 
domain  can  result  in  improved  responses  (Fig.  8  (b.ii))  reducing  TEDU 2  to  1040  Jm-2.  The  sponge 
layer  thus  allows  for  significant  leeway  in  boundary  positioning  but  it  does  not  necessarily  entirely 
eliminate  issues  in  the  response.  Note  how  in  the  two  open  boundaries  domain  the  amphidrome  in 
the  Southern  Ocean  below  Australia  is  located  near  where  the  third  open  boundary  is.  The  sponge 
acts  to  push  this  amphidrome  away  from  the  boundary  because  the  three  open  boundaries  IndWPac 
model  solution  and  the  TPX08  solution  are  incompatible  here.  This  indicates  reliance  on  internal 
dissipative  mechanisms  to  ensure  compatibility  with  each  other. 

5.1.3  Discussion 

What  to  make  of  the  dramatic  results  to  the  solution  due  to  boundary  placement  and  condition 
types?  Firstly,  even  though  there  should  not  be  an  amphidrome  right  next  to  the  boundary  ac¬ 
cording  to  the  TPX08  solution  (although  the  elevations  are  still  fairly  small),  our  model  without 
adequate  internal  dissipative  effects  expects  there  to  be  one.  Instabilities  and  problems  may  arise 
near  amphidromes  because  a  physically  incorrect  solution  that  satisfies  the  governing  equations  can 
be  obtained.  Mathematically,  both  boundary  conditions  and  initial  conditions  are  required  to  get 
the  correct  solution.  Instead,  the  method  commonly  adopted  (including  this  study)  is  to  impose 
the  elevation  boundary  conditions  and  ramp  up  the  system  from  a  completely  zero  state.  We  found 
from  a  simple  test  case  that  when  internal  dissipative  effects  are  low,  ramping  generates  spurious 
modes  that  persist  for  very  long  periods  of  time.  Secondly,  clearly  the  Indian  Ocean  and  the  Aus¬ 
tralian/Indonesian  marginal  seas  appear  very  sensitive  to  fluctuations  in  fluxes  on  the  boundary. 
As  noted  by  Zahel  and  Muller  (2005)  in  their  study  of  free  barotropic  oscillations,  the  11.65  hour 
resonant  mode  shows  a  very  similar  pattern  to  the  lunar  semi-diurnal  tides  in  the  Indian  Ocean. 
Furthermore,  the  energy  density  of  this  mode  is  2.2  times  the  global  average  in  the  Indian  Ocean. 
Hence  the  resonant  nature  of  the  lunar  semi-diurnal  tides  in  this  basin  causes  TED  to  increase 
wildly  in  response  to  the  poor  boundary  conditions. 

However,  one  of  the  main  issues  with  using  the  absorption-generation  sponge  layer  is  the  reliance 
on  the  reference  solution.  In  particular,  the  fluxes  obtained  from  TPXO8  are  less  likely  to  be 
as  accurate  or  compatible  with  the  IndWPac  model  as  the  TPX08  elevations  (because  TPX08 
only  assimilates  elevations).  We  found  an  interesting  case  in  our  model  where  the  solution  near 
the  southwest  boundary  degrades  qualitatively  versus  the  TPX08  solution  when  implementing  the 
sponge  layer  for  the  current  two  boundary  domain  (Fig.  8  (a.ii)).  It  is  also  shown  in  §5.3  that 
dissipation  in  the  IndWPac  model  is  slightly  negative  over  a  wide  region  of  the  Southern  Ocean  which 
could  be  due  to  the  boundary  conditions.  Perhaps  differences  in  fluxes  are  accentuated  due  to  the 
skewness  of  the  projection  at  high-latitudes.  Additionally,  in  the  Southern  Ocean  progressive  waves 
are  allowed  to  freely  propagate  around  the  globe  suggesting  that  fluxes  exiting  the  boundaries  here 
would  otherwise  be  allowed  to  enter  the  opposite  boundary  eventually  leading  to  some  equilibrium 
condition  (that  may  be  somewhat  different  than  TPX08)  in  a  global  model. 
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5.2  Bathymetry 

Bathymetry  is  a  boundary  condition  for  oceanic  models  hence  its  importance  to  the  solution  is 
clear.  Recent  years  have  shown  marked  improvements  in  global  bathymetric  databases  such  as 
SRTM15JPLUS  and  GEBCO-2014.  This  section  begins  by  outlining  the  effects  of  using  one  of  these 
databases  over  the  other,  followed  by  effects  between  SRTM15JPLUS  and  our  more  comprehensive 
bathymetric  data  (Table  1).  The  section  concludes  with  a  discussion  on  the  results  and  implications. 

5.2.1  Comparisons  between  Global  Bathymetric  Databases 

Current  global  bathymetric  databases  SRTM15.PLUS  and  GEBCO-2014  are  good  enough  that 
they  do  allow  us  to  obtain  mean  RMS  discrepancies  ~3  cm  for  the  M2  tidal  wave  in  the  deep 
ocean.  Nevertheless,  there  is  still  a  reasonable  level  of  uncertainty  between  them  (Fig  9).  For 
example,  on  the  abyssal  hills  the  use  of  statistical  roughness  (Goff  and  Arbic,  2010)  to  calculate  a 
new  bathymetry  set  (Timko  et  ah,  2017)  has  been  undertaken  to  account  for  the  effective  coarseness 
of  satellite  altimetry  derived  bathymetry  (note  that  the  SRTM15-PLUS  database  used  here  is  a 
combination  with  SRTM30-PLUS  that  contains  the  synthetic  realization  of  the  abyssal  hill  roughness 
but  GEBCO_2014  is  without  it).  To  investigate  the  effects  of  the  global  bathymetric  databases  we 
compute  one  simulation  using  SRTM15-PLUS  with  abyssal  hill  roughness  everywhere  ( SRTM  IT 
-f  SC)  and  another  using  GEBCO-2014  everywhere  (GEBCO  IT  +  SC).  RMS  differences  between 
the  simulations  for  the  M2  tidal  wave  are  plotted  in  Fig  9.  Optimal  internal  tide  dissipation  factors 
(IT)  and  spatially  constant  Cj  =  2.5xl0-3  (SC)  are  employed  for  both. 

Along  the  ocean  ridges  that  include  a  synthetic  realization  of  the  abyssal  hill  roughness,  the 
normalized  bathymetric  differences  range  between  5-25%  except  for  the  Southwest  Indian  Ocean 
Ridge  where  it  can  differ  by  more  than  50%  in  spots.  Despite  this,  M2  RMS  differences  in  deep  water 
do  not  exceed  2  cm  anywhere  except  cast  of  Australia  and  New  Guinea,  i.e.  the  RMS  differences 
between  responses  resulting  from  the  two  global  bathymetric  databases  tend  to  be  much  less  than 
RMS  discrepancies  between  Comp  +  IT  +  SV  model  setup  and  TPX08  (Fig  4  (c.i)).  Note  that 
although  differences  in  the  bathymetry  should  change  the  internal  tide  dissipation  matrix,  we  used 
the  same  dissipation  matrix  as  the  Comp  -h  IT  +  SV  model  setup  for  both  simulations  to  help 
identify  strictly  bathymetric  effects. 

Major  normalized  bathymetric  and  RMS  differences  are  unsurprisingly  found  in  shallow  waters 
such  as  the  South  China  Sea,  Bass  Strait,  Yellow  Sea,  Sea  of  Okhotsk,  Ganges  Delta,  northern 
Andaman  Sea,  Persian  Gulf,  and  the  Gulfs  of  Khambhat  and  Kutch.  The  most  astounding  RMS 
differences  are  located  on  the  northern  Australian  shelves,  the  Taiwan  Strait  and  the  Seto  Inland 
Sea.  In  particular,  the  RMS  differences  in  the  Coral  Sea  and  around  Torres  Strait  are  larger  than 
the  discrepancies  between  the  Comp  -h  IT  -h  SV  model  setup  and  TPX08.  According  to  the  sources 
of  GEBCO-2014,  the  2009  Australian  Bathymetry  and  Topography  Grid  (Whiteway,  2009)  is  used 
around  the  Australian  continent.  Within  this  dataset,  some  of  the  nearshore  bathymetry  is  made 
up  of  multibeam,  Laser  Airborne  Depth  Sounder,  and  nautical  charts,  with  the  rest  based  on  1  arc 
min  and  2  arc  min  ETOPO  satellite  derived  bathymetry.  In  comparison,  SRTM  15-PLUS  is  said 
to  include  50  m  multibeam  datasets  from  2012  as  well  as  the  Deepreef  Explorer  GBR  dataset  from 
2010  (both  newer  than  the  2009  Australian  Bathymetry  and  Topography  Grid). 

5.2.2  Comparisons  between  SRTM15-PLUS  and  Local  High-Resolution  Bathymetry 

Another  issue  with  a  global  bathymetric  dataset  such  as  SRTM  15-PLUS  is  that  on  the  shelf  and 
nearshore  the  resolution  can  be  too  coarse  and  it  may  contain  holes  in  the  bathymetry,  thus  it 
is  not  completely  reliable  for  accurate  regional  simulations.  This  section  outlines  the  differences 
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Figure  9:  Differences  between  SRTM15-PLUS  with  abyssal  hill  roughness  (Goff  and  Arbic,  2010)  and 
GEBCO.2014  global  bathymetric  databases,  (a)  Normalized  differences  in  the  bathymetry,  (b)  M2  RMS 
differences  in  the  responses  between  the  GEBCO  IT  +  5(7,  and  SRTM  IT  +  SC  model  setups  (see  Table  4 
for  description  of  model  setups). 


between  SRTM15-PLUS  and  our  more  comprehensive  bathymetric  data  (Table  1)  containing  local 
high-resolution  datasets.  We  focus  on  three  regions:  the  East  China  Sea  including  the  Yellow  Sea 
and  southern  Japan;  the  South  China  Sea  including  the  Philippines  Seas  and  north  Java  Sea;  and 
the  Coral  Sea  including  the  Torres  Strait  (Fig  10). 

All  nearshore  areas  in  the  East  China  Sea  show  significant  normalized  bathymetric  differences 
aside  from  Hong  Kong  which  contains  our  smallest  element  sizes.  Due  to  numerous  spurious  large 
depths  in  the  SRTM15_PLUS  dataset  in  this  region  we  had  to  replace  this  area  with  the  local  high- 
resolution  bathymetry  in  order  to  avoid  instabilities  due  to  violation  of  the  CFL  condition.  This  is 
not  thought  to  have  a  large  effect  on  the  results  of  the  comparisons  between  the  bathymetric  datasets 
and  the  conclusions  that  we  draw  from  them.  The  simulations  show  large  RMS  differences  of  the  M2 
tidal  wave  in  the  Taiwan  Strait  and  Seto  Inland  Sea  but  interestingly  they  are  not  as  large  as  those 
between  GEBCO  and  SRTM  15-PLUS,  nor  are  differences  in  the  Gulf  of  Tonkin  as  pronounced.  It 
should  be  mentioned  that  we  did  notice  reduced  discrepancies  at  tide  gauges  in  the  Seto  Inland 
Sea  when  using  the  high-resolution  bathymetry.  It  is  a  complicated  region  with  many  small  islands 
and  channels  and  requires  accurate  connectivity  of  the  energy  fluxes  to  improve  results.  Most  of  the 
RMS  differences  in  the  Yellow  Sea  between  SRTM15JPLUS  and  the  local  high- resolution  bathymetry 
are  larger  than  those  between  GEBCO  and  SRTM15-PLUS,  in  particular  north  of  Shanghai  and  in 
the  Incheon  area.  However,  generally  these  differences  are  smaller  than  the  discrepancies  between 
TPX08  and  the  Comp  +  IT  +  SV  model  setup. 

Despite  widespread  normalized  bathymetric  differences  in  the  Philippines  and  the  region  between 
the  South  China  Sea  and  Java  Sea,  fairly  small  RMS  differences  in  the  M2  tidal  wave  result  aside 
from  a  few  channels  near  Singapore,  and  in  the  gulfs  near  Batang  Lupar  and  North  Kalimantan, 
both  on  Borneo.  This  is  perhaps  partly  because  the  M2  amplitudes  are  relatively  small  in  the  region 
between  the  South  China  Sea  and  Java  Sea.  The  larger  Ki  (not  shown)  shows  up  to  10-30  cm 
RMS  difference  in  the  area  south  of  Singapore  but  is  not  notable  elsewhere.  In  the  two  gulfs  on 
Borneo,  which  have  fairly  large  M2  tidal  ranges  (up  to  1.7  m  in  the  gulf  near  Batang  Lupar),  the 
differences  result  not  only  from  local  high-resolution  bathymetric  datasets  but  due  to  hand-edits  vis- 
a-vis  FUGAWI  navigational  charts.  Areas  like  the  gulf  near  Batang  Lupar  can  be  extremely  sensitive 
to  bathymetry  particularly  deep  in  the  gulf  where  the  v-shape  concentrates  the  tidal  energy.  The 
effect  of  our  hand-edits  was  to  deepen  the  area  near  Batang  Lupar  allowing  the  tidal  range  to  reach 
close  to  the  measured  one. 

The  final  region  is  the  Coral  Sea  which  demonstrably  shows  large  widespread  normalized  bathy¬ 
metric  differences.  This  is  slightly  perplexing  as  SRTM1 5-PLUS  should  include  the  Deepreef  Ex¬ 
plorer  GBR  dataset  according  to  their  references  (ftp://topex.ucsd.edu/pub/srtml5_plus/). 
However,  the  SRTM1 5-PLUS  version  used  in  this  study  does  not  seem  to  have  it  incorporated. 
The  bathymetry  in  the  Torres  Strait  and  Papua  New  Guinea  region  is  also  very  different  from 
SRTM1 5_PLUS,  where  we  have  actually  used  GEBCO-2014  in  our  more  comprehensive  bathymetric 
dataset.  This  is  because  GEBCO-2014  matches  rather  well  with  Deepreef  Explorer  GBR  at  the 
interface  whereas  SRTM  15-PLUS  does  not.  The  resulting  M2  RMS  differences  shown  here  are  cer¬ 
tainly  large  and  often  well  exceed  those  discrepancies  between  our  model  and  TPX08.  In  fact,  in 
most  of  the  Coral  Sea  the  Comp  -h  IT  -h  SV  model  setup  is  performing  rather  well  with  respect 
to  TPX08  which  could  be  largely  attributed  to  the  Deepreef  Explorer  GBR  dataset.  Fairly  large 
discrepancies  for  our  model  against  TPX08  and  tide  gauges  are  still  present  near  the  Torres  Strait 
where  opposing  M2  energy  fluxes  meet  over  the  strait  (Fig  11,  the  definition  of  energy  flux  density, 
P  is  presented  in  §5.3).  The  residual  discrepancy  here  is  likely  a  combination  of  the  remaining  un¬ 
certainties  in  the  GEBCO-2014  bathymetry  (based  on  2009  Australian  Bathymetry  and  Topography 
Grid),  and  bottom  friction  dissipation. 
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Figure  10:  Differences  between  SRTM15.PLUS  with  abyssal  hill  roughness  (Goff  and  Arbic,  2010)  and 
our  comprehensive  bathymetric  data  (Table  1)  containing  local  high-resolution  datasets,  (i)  Normalized 
differences  in  the  bathymetry,  (ii)  M2  RMS  differences  in  the  responses  between  SRTM  +  IT  +  SC  and 
Comp  +  IT  +  SC  model  setups  (see  Table  4  for  description  of  model  setups),  (a)  East  China  Sea,  (b)  South 
China  Sea,  (c)  Coral  Sea. 
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Figure  11:  Energy  flux  density  P  of  the  M2  tidal  wave  over  the  Torres  Strait. 


5.2.3  Discussion 

The  effect  of  different  bathymetric  datasets  is  not  shown  to  be  an  important  factor  in  the  deep 
ocean,  but  on  the  shelf  and  nearshore  there  are  certain  regions  where  bathymetry  plays  a  large  role. 
This  was  also  highlighted  in  terms  of  the  global  RMS  discrepancies  presented  in  §4.5  (summarized 
in  Table  4).  In  some  cases  the  RMS  differences  between  simulations  using  different  bathymetric 
datasets  were  greater  than  discrepancies  between  TPX08  and  the  Comp  +  IT  +  SV  model  setup, 
particularly  between  the  two  global  bathymetric  datasets.  Our  more  comprehensive  bathymetric 
data  and  SRTM  15-PLUS  are  mostly  the  same  except  nearshore  in  certain  regions  which  was  likely 
the  reason  for  smaller  RMS  differences  in  general. 

As  discussed  in  Zaron  (2017),  bathymetry  has  the  potential  to  control  the  tidal  elevation  partic¬ 
ularly  through  resonant  effects  that  are  in  general  nonlocal.  For  example,  the  effect  of  the  Deepreef 
Explorer  GBR  high-resolution  bathymetry  was  to  change  the  M2  elevation  over  a  large  area  beyond 
Athe  Coral  Sea  out  into  the  deep  ocean.  Conversely  only  small  changes  were  noted  throughout  the 
region  between  the  South  China  Sea  and  Java  Sea.  More  locally,  in  a  resonant  basin  such  as  the  gulf 
near  Batang  Lupar,  hand-edits  of  the  bathymetry  based  on  FUGAWI  navigational  charts  allowed 
the  tidal  elevation  to  reach  close  to  the  measured  M2  amplitude.  Clearly,  greater  availability  and 
quality  of  nearshore  and  shelf  bathymetry  has  the  potential  to  greatly  improve  the  modeling  not 
only  of  tides  but  all  shallow  water  flows. 

With  regards  to  the  effect  of  the  bathymetry  in  the  deep  ocean,  it  should  be  noted  that  bathymetry 
will  affect  internal  tide  dissipation  as  the  dissipation  matrix  is  based  on  topographic  depths  and 
slopes.  So  technically,  the  deep  ocean  may  be  more  impacted  than  is  shown  here  taking  this  aspect 
into  account.  Nevertheless,  addition  of  the  abyssal  hill  roughness,  for  example,  has  only  marginally 
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increased  the  fidelity  of  ocean  models  (Buijsman  et  al.,  2015;  Timko  et  al.,  2017).  In  3D  baroclinic 
models  (Timko  et  al.,  2017)  this  can  be  explained  in  part  by  limitations  of  resolution.  On  the  other 
hand,  2D  barotropic  models  such  as  IndWPac  can  achieve  high-resolution  over  a  wide-scale,  but 
may  be  somewhat  limited  by  the  underlying  assumptions  of  internal  tide  dissipation  no  matter  the 
bathymetric  data.  Greater  discussion  on  this  aspect  is  presented  in  the  next  section. 

5.3  Internal  Tide  Dissipation 

In  the  two  internal  tide  dissipation  parameterizations  (§3.3),  it  is  necessary  to  calibrate  a  global 
amplification  factor  due  to  unknowns  involved  with  the  resolution  of  the  bathymetric  data  and  the 
way  in  which  dissipation  that  is  overestimated  at  supercritical  slopes  is  handled,  due  to  their  linear 
assumptions.  In  addition  to  finding  these  amplification  factors,  this  section  discusses  the  differences 
between  the  two  parameterization  methods,  introduces  multiplier  coefficients  due  to  semi-diurnal 
resonance  in  the  Luzon  Strait,  and  concludes  with  some  final  remarks  on  reasons  for  differences 
between  the  methods  and  remaining  issues  for  the  parameterization  of  internal  tide  dissipation. 


Dissipation  [GW] 

Figure  12:  Total  dissipation  TD  versus  RMS  discrepancy  against  TPX08  Dtpx  in  deep  water  ( h  >  500  m) 
for  different  internal  tide  parameterization  methods  ( Nycander  and  Directional)  and  amplification  factors 
( CnVc  and  Cdw  respectively)  with  2nd  order  polynomial  fits;  (a)  M2  tidal  wave,  (b)  Ki  tidal  wave.  Nycander 
method  without  positive  definite  (PD)  or  Nb  spatial  averaging  ( Nbav )  modifications,  and  TPXO8  dissipation 
in  deep  water  are  shown  for  reference. 
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5.3.1  Calibrating  Amplification  Factors 

We  begin  by  trying  to  determine  the  optimal  values  of  C^yc  and  CW  for  the  IndWPac  model  before 
comparing  the  performance  of  both  parameterization  methods.  This  is  evaluated  by  looking  at  Dtpx 
for  the  M2  and  Ki  tidal  waves  in  deep  water  ( h  >  500  m)  with  a  spatially  constant  bottom  friction, 
Cf  =  2.5xl0~3.  To  a  lesser  degree  we  are  also  interested  in  the  total  dissipation  TD  of  individual 
tidal  waves: 

TD  =  JJ  (W  -  V  P)  dA  (16) 

w  =  9JtL  [  u' V(7?£Q  +  vs al)  dt  (17) 

p=9£oh£ur,dt  (18) 

where,  W  is  the  work  rate,  and  P  is  the  energy  flux.  Since  the  numerical  calculation  of  V  •  P 
with  finite  precision  is  very  noisy,  in  this  work  the  area-integral  is  computed  using  the  divergence 
theorem.  Comparisons  of  Dtpx  versus  TD  in  deep  water  for  the  individual  tidal  waves  are  shown  in 
Fig.  12  using  four  different  values  of  amplification  factors  for  each  method. 

Regarding  the  Ki  tidal  wave,  amplification  factors  slightly  smaller  than  the  values  tested  appear 
optimal  in  both  parameterization  methods.  However,  we  focused  on  the  results  of  the  M2  tidal  wave 
to  determine  the  optimal  amplification  factors.  This  led  to  C^yc  &  2.9  and  Cjr>jr  «  0.22  based  on  a 
second-order  polynomial  best  fit.  For  comparison,  Buijsman  et  al.  (2015)  determined  C^yc  ~  2.75, 
which  is  in  fairly  close  agreement.  Furthermore,  two  modifications  were  applied  to  the  Nycander 
method.  The  first  modification  simply  ensured  that  the  dissipation  matrix  was  positive  definite 
(PD).  A  second  modification  involved  applying  Gaussian  smoothing  of  Nb  (to  obtain  a  variable 
denoted  as  Nbav)  using  the  same  radius  and  scaling  as  the  convolution  integral  for  J,  the  argument 
being  that  nonlocal  effects  of  buoyancy  may  be  just  as  important  as  nonlocal  effects  of  topography. 
Both  modifications  create  more  dissipation  and  makes  Dtpx  for  M2  smaller  for  the  same  value  of 
pNyc  (Fig.  12),  with  the  PD  modification  having  the  largest  effect. 

5.3.2  Differences  between  Parameterization  Methods 

Two  opposing  outcomes  result  from  the  comparison  between  the  Nycander  and  Directional  methods. 
For  both  constituents  the  Nycander  method  leads  to  slightly  smaller  values  of  Dtpx ,  while  TD  at 
the  optimal  amplification  factor  matches  TPX08  in  deep  water  more  closely  for  the  Directional 
method.  For  the  reason  that  TPX08  can  be  reliably  validated  for  elevations  but  not  for  dissipation 
we  are  inclined  to  the  prefer  the  Nycander  method,  thus  it  is  incorporated  into  our  best  model 
results  presented  in  §4.  It  is  however  worth  pointing  out  that  the  difference  between  the  methods  is 
no  more  that  15  mm  in  Dtpx  for  M2,  thus  the  Directional  method  can  be  considered  a  very  useful 
parameterization  in  its  own  right  -  not  least  because  it  can  be  quickly  calculated  and  introduced  to 
a  numerical  model. 

With  regards  to  total  dissipation,  the  optimal  Nycander  method  results  in  24%  greater  M2  TD 
compared  with  TPX08,  which  is  in  very  close  agreement  to  the  HYCOM  global  model  (23%  larger 
than  TPX08)  (Buijsman  et  al.,  2015).  Buijsman  et  al.  (2015)  notes  that  the  TPX08  dissipation 
rates  are  diffused  over  large  areas  in  comparison  to  the  parameterized  internal  tide  dissipation  in 
their  model.  In  that  sense  it  is  somewhat  unclear  how  reliable  TPX08  dissipation  in  deep  water  may 
be.  Issues  with  tidal  dissipation  in  global  data-assimilated  models  have  been  previously  highlighted 
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(Lyard  et  al.,  2006;  Le  Provost  and  Lyard,  1997),  and  there  are  large  regions  of  negative  dissipation 
rates  when  computing  this  with  the  TPX08  solutions. 

To  get  a  sense  of  the  global  differences  in  the  way  the  two  methods  dissipate  energy  and  affect 
the  response,  plots  of  the  dissipation  density  (W  —  V  •  P)  for  M2  between  the  optimal  Nycander 
and  Directional  methods  are  shown  in  Fig.  13.  The  spatial  distribution  is  similar  to  previous  studies 
(e.g  Green  and  Nycander,  2013;  Buijsman  et  al.,  2015).  In  the  deeper  regions  such  as  the  abyssal 
hills  (namely  the  West  Indian  Ocean  Ridge  and  Western  Pacific  Trench),  the  general  observation 
is  that  the  Nycander  method  focuses  dissipation  towards  the  center  peaks  of  the  ridges  whilst  the 
Directional  method  tends  to  spread  dissipation  over  the  width  of  the  ridge.  This  makes  sense  because 
the  Directional  method  is  a  function  of  the  square  of  the  local  topographic  gradient  meaning  that 
dissipation  may  be  triggered  locally  more  easily  than  the  Nycander  method  which  requires  both  a 
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Figure  13:  Dissipation  density  of  the  M2  tidal  wave  comparing  internal  tide  dissipation  methods;  (a) 
Nycander  method  (PD  and  Nbavg  corrections)  with  CnVc  —  2.90,  (b)  Directional  method  with  CDir  —  0.22. 
hatched  regions  indicate  sponge  layer. 
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large  local  topographic  gradient  and  a  large  gradient  of  J  (which  will  be  largest  in  the  center  or 
ridges  as  it  is  a  measure  of  the  overall  topographical  change  in  the  vicinity).  This  is  true  too  for 
island  chains  where  the  Nycander  method  tends  to  give  greater  dissipation  here  because  even  if  the 
local  topographic  gradients  are  not  large,  the  topographical  change  in  the  general  vicinity  of  the 
islands  might  be. 

Depth-wise  the  characteristics  of  local  and  large-scale  topographic  gradients  tend  to  translate 
into  the  Nycander  method  creating  more  dissipation  in  shallower  depths  (Fig.  14).  Both  methods 
give  large  amounts  of  dissipation  in  the  3000  -  4000  m  range  corresponding  to  abyssal  hills,  but 
the  Directional  method  dissipation  drops  off  more  rapidly  with  decreasing  depth  to  match  TPX08 
surprisingly  well  in  shallow  depths  <  500  m.  As  mentioned  above,  it  remains  unclear  whether  the 
overestimate  in  dissipation  from  both  schemes  in  the  500  -  4000  m  depth  range  (and  underestimate 
in  the  4500  -  6000  m  range)  is  a  major  cause  of  concern  or  that  it  simply  reflects  the  coarseness 
of  the  data-assimilated  model.  Buijsman  et  al.  (2015)  found  similar  trends  (to  this  study)  for  their 
global  model. 

The  differences  in  amplitudes  of  the  M2  tidal  wave  between  the  Nycander  and  Directional  meth¬ 
ods  are  illustrated  in  Fig.  15.  There  is  a  clear  divide  between  amplitudes  in  the  Indian  Ocean  and 
those  in  the  Western  Pacific  Ocean.  This  indicates  disparity  in  the  way  the  M2  tides  are  balanced 
between  basins  depending  on  the  method.  The  Nycander  method  dissipates  more  in  shallower  depths 


Total  Dissipation  [GW] 

Figure  14:  Bathymetric  depth  versus  total  dissipation  TD  (summed  in  100  m  depth  bins)  for  the  M2  tidal 
wave  comparing  TPXO8  and  the  internal  tide  dissipation  parameterization  methods  (modified  Nycander  and 
Directional)  with  their  optimal  amplification  factors.  Nycander  method  without  the  Nb  spatial  averaging 
modification  ( Nbav )  is  shown  for  reference. 
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(Fig.  14),  so  it  is  perhaps  unsurprising  that  the  amplitudes  will  be  smaller  in  the  Western  Pacific 
basin  which  contains  many  shallow  shelves  and  island  chains.  Moreover,  the  energy  flux  density 
P  of  the  M2  tidal  wave  predominantly  flows  from  the  Indian  Ocean  into  the  Western  Pacific  basin 
through  the  Indonesian  Seas,  and  to  a  lesser  extent  through  the  Malacca  Strait  (Fig.  16).  Due  to  the 
greater  dissipation  in  shallow  depths  as  it  passes  through  the  island  chains  and  shallow  seas,  more 
energy  remains  on  the  Indian  Ocean  side  instead  of  flowing  into  the  Western  Pacific  basin  compared 
with  the  Directional  method. 

5.3.3  Multiplier  due  to  Semi-diurnal  Resonance  over  the  Luzon  Strait 

A  local  improvement  to  internal  tide  dissipation  was  applied  over  the  Luzon  Strait.  The  Luzon 
Strait  controls  the  amount  of  energy  into  the  South  China  Sea  for  both  the  M2  and  Ki  tidal  waves 
(Fig.  16).  Here,  the  amplitude  of  Ki  is  mostly  larger  than  M2  (Fig  4).  Because  the  separation  of 
the  double-ridge  topography  in  the  strait  is  similar  to  the  semi-diurnal  internal  tide  wavelength,  it 
has  been  shown  that  resonance  dramatically  increases  the  barotropic  to  baroclinic  energy  conversion 
(Buijsman  et  al.,  2014).  In  comparison  to  the  sum  of  the  two  ridges  considered  separately,  the 
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Figure  15:  M2  tidal  wave  amplitude  differences  between  the  two  internal  tide  dissipation  methods  (modified 
Nycander  method  with  CnVc  =  2.9  minus  the  Directional  method  with  Coir  =  0.22). 
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Figure  16:  Energy  flux  density  P  of  the  (a)  M2  and  (b)  Ki  tidal  waves  in  the  marginal  seas  separating  the 
Indian  Ocean  from  the  Western  Pacific  Ocean. 
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double-ridge  produced  up  to  as  much  as  four  times  the  energy  conversion  for  the  first-internal  mode 
(Buijsman  et  al.,  2014).  Since  our  internal  tide  dissipation  parameterizations  do  not  include  such 
resonance  behavior  we  deemed  it  appropriate  to  apply  a  multiplier  to  the  amplification  coefficients 
in  the  region  19.5°-21.5°N  and  120°-122.5°E.  The  multiplier  coefficients  C Luzon  were  defined  using 
a  skewed  Gaussian  curve  as  a  function  of  latitude  <j>  to  approximate  the  data-points  presented  in 
Buijsman  et  al.  (2014): 


with  cll  ~  5.0,  crL  =  0.3,  oll  =  -1.0,  and  £l  =  (j)  —  20.9°N.  Cluzou  reaches  a  maximum  of  4.24  at 
20.75°N. 

Fig.  17  shows  the  amplitude  differences  for  the  M2  and  Kj  tidal  waves  in  the  South  China  Sea 
and  surrounds  between  the  optimal  Nycander  method  with  and  without  the  multiplier  coefficients 
C Luzon-  As  expected,  the  increase  to  the  internal  tide  coefficients  reduces  the  amplitudes  inside 
most  of  the  South  China  Sea  for  both  constituents.  The  decrease  is  on  the  order  of  0.5-1  cm  for 
M2  in  most  areas  with  a  few  pockets  of  2.5-5  cm  reductions.  Additionally,  the  blockage  increases 
amplitudes  slightly  to  the  east  of  the  strait  and  in  the  Sulu  Sea.  The  Ki  amplitude  is  decreased 
in  the  South  China  Sea  by  1-2.5  cm  almost  uniformly.  Furthermore,  due  to  the  blockage  at  the 
Luzon  Strait  more  of  energy  flux  is  now  diverted  down  into  the  Indonesian  seas  (Fig.  16)  increasing 
amplitudes  uniformly  by  0.5-1  cm. 

M2  RMS  discrepancies  at  coastal  tide  gauges  are  generally  decreased  on  the  order  of  1  cm 
within  the  South  China  Sea  due  to  the  Cluzou  multiplier  coefficients.  The  discrepancy  is  increased 
slightly  in  the  Sulu  Sea  and  near  the  Taiwan  Strait.  The  changes  in  Ki  discrepancies  do  not 
follow  a  clear  pattern  aside  from  in  the  Gulf  of  Thailand  and  Celebes  Sea  regions  (decrease  and 
increase  respectively).  The  RMS  discrepancies  against  tide  gauges  and  TPX08  in  the  plotted  region 
(Fig.  17)  are  summarized  in  Table  5.  Overall  only  small  decreases  in  discrepancies  are  found  when 
using  the  Cluzou  multiplier  coefficients.  However,  the  discrepancies  compare  rather  favorably  to  a 
local  hydrodynamic  model  for  just  the  South  China  Sea  region  (Green  and  David,  2013)  (we  get  8.5 
cm  and  5.0  cm  RMSE  for  M2  and  Ki  versus  TPX08  respectively,  Green  and  David  (2013)  quotes 
9  cm  and  10  cm).  In  that  study  Cf  had  to  be  raised  to  an  unphysical  value  of  0.01  to  achieve 


Table  5:  The  mean  RMS  discrepancies  (units:  cm)  versus  coastal  tide  gauges  Dtg  and  TPX08  Dtpx  (in 
all  depths)  of  the  M2,  Ki,  and  the  total  free  surface  (up  to  all  eight  major  constituents)  within  each  of 
the  regions  plotted  in  Figs.  17  and  19.  Standard  deviations  in  parentheses.  See  Table  4  for  model  setup 
descriptions. 


Region 

Model 

M2  nKl  all  _^M2  r>Kl  all 

UtQ  JJtQ  UtQ  Ulvx  utvx  utvx 

scs 

Comp  +  IT(LZ)  +  SV 
Comp  +  IT(NoLZ)  +  SV 

7.35  (6.75)  6.92  (5.88)  14.3  (10.2)  3.76  (7.56)  3.63  (3.37)  6.70  (10.4) 

7.55  (6.63)  6.96  (5.73)  14.2  (10.0)  4.07  (7.53)  3.74  (3.56)  6.90  (10.4) 

YS 

Corny  +  IT  +  SC 

Comp  +  IT  +  SV 

17.2  (17.7)  3.74  (3.74)  20.8  (20.3)  8.58  (16.9)  1.70  (9.53)  10.3  (18.4) 

12.3  (15.6)  3.04  (3.31)  16.1  (17.7)  6.51  (17.1)  1.41  (9.81)  8.43  (18.2) 

JS 

Comp  -h  IT  -f-  SC 

Comp  +  IT  +  SV 

10.5  (9.50)  5.98  (5.41)  16.5  (11.6)  6.02  (7.11)  5.47  (7.08)  10.6  (7.39) 

11.5  (10.8)  6.54  (5.39)  17.9  (12.6)  5.97  (7.05)  6.60  (6.66)  11.5  (7.58) 

TAS 

Comp  ~h  IT  -h  SC 

Comp  +  IT  +  SV 

18.8  (24.6)  9.49  (8.16)  25.9  (28.2)  6.46  (9.68)  4.03  (8.50)  9.71  (11.1) 

19.9  (24.8)  10.2  (8.63)  27.6  (28.6)  6.81  (9.86)  4.28  (8.62)  10.3  (11.7) 

*SCS:  South  China  Sea  region  plotted  in  Fig.  17.  LZ  refers  to  the  use  of  multiplier  coefficients,  Cluzou  from  (19), 

over  the  Luzon  Strait.  NoLZ  is  without  applying  Cluzou 

*YS:  Yellow  Sea  and  southern  Japan  region  plotted  in  Fig.  19  (i) 

*JS:  Area  between  the  Java  Sea  and  South  China  Sea  plotted  in  Fig.  19  (ii) 

*TAS:  Timor  and  Arafura  Seas  region  plotted  in  Fig.  19  (iii) 


34 


25* 

20* 

15- 

ID' 

5’ 

25‘ 

20* 

15' 

IQ¬ 
S’ 

100'  105-  110*  115’  120’  125* 

Figure  17:  (a)  M2  and  (b)  Ki  amplitude  differences  in  the  South  China  Sea  and  surrounds  when  using 
the  multiplier  coefficients  C Luzon  over  the  Luzon  Strait  (case  with  C Luzon  from  (19)  minus  case  without 
C Luzon  applied).  Circles  indicate  the  change  in  RMS  discrepancies  at  coastal  tide  gauges  (negative  indicates 
reduction  in  discrepancy  when  using  C Luzon)- 


optimal  results  for  M2.  It  was  speculated  that  this  is  because  the  internal  tide  dissipation  rates 
in,  e.g.  the  Luzon  Strait,  are  underestimated  for  M2.  This  may  have  been  the  case  in  Green  and 
David  (2013),  but  the  results  shown  here  indicate  that  using  C Luzon  multiplier  coefficients  due  to 
semi-diurnal  resonance  from  the  double-ridge  only  provides  a  marginal  positive  effect.  One  reason 
for  the  disparity  between  the  two  models  could  be  that  fluxes  to  the  east  of  the  Luzon  Strait  may 
be  quite  different  since  the  Green  and  David  (2013)  model  is  driven  by  the  TPX08  solution  at 
lateral  boundaries  close  to  the  strait  whereas  our  boundaries  are  much  further  away.  Differences 
in  bathymetry  (SRTM15.PLUS  and  the  inclusion  of  local  high-resolution  bathymetry  in  IndWPac, 
versus  GEBCO.2014  in  Green  and  David  (2013))  also  likely  play  a  role  in  improved  responses  for 
the  IndWPac  model. 
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5.3.4  Discussion 


As  a  result  of  the  dissipation  dynamics  we  found  that  the  Directional  method  gave  better  results  in 
the  Western  Indian  Ocean,  but  elsewhere  the  Nycander  method  was  generally  preferable.  It  is  worth 
noting  that  the  M2  amplitude  differences  (Fig  15)  in  the  Indian  Ocean  between  the  two  methods 
ranges  between  0.5-2. 5  cm,  and  more  in  the  Bay  of  Bengal.  This  is  a  rather  large  amplitude  difference 
in  the  deep  ocean  since  Dtpx  =  2.9  cm,  indicating  that  there  is  some  scope  to  improve  deep  water 
solutions  further  through  better  parameterization  of  internal  tide  dissipation.  Perhaps  some  of  the 
remaining  issues  for  the  Nycander  method  can  be  explained  by  the  fact  that  bathymetry  is  still  rather 
uncertain  and  coarse  in  much  of  the  deep  ocean  (most  of  the  Indian  Ocean  is  deep  with  very  narrow 
shelves,  and  internal  tide  generation  is  an  important  dissipation  mechanism).  Furthermore,  internal 
tide  dissipation  in  shallow  regions  is  less  reliable  because  of  larger  tidal  velocities  and  uncertainties, 
and  there  is  a  greater  chance  of  the  flow  being  supercritical  (Melet  et  al.,  2013).  In  fact,  one  of  the 
main  effects  of  the  Nbav  modification  to  the  Nycander  method  was  to  move  some  dissipation  away 
from  shallow  regions  into  deeper  regions  (Fig.  14).  More  investigation  into  the  parameterization  of 
internal  tide  generation  in  shallow  regions  is  warranted  especially  because  the  shallow  areas  of  the 
Indonesian  seas  provide  a  critical  connection  between  the  basins,  a  region  that  has  created  issues 
previously  (Melet  et  al.,  2013). 

Finally,  it  is  worth  highlighting  that  the  internal  tide  dissipation  matrices  used  here  have  been 
derived  based  on  the  M2  tidal  wave.  That  is,  lj  =  2tt/12.42  rad/hr  is  used  in  (4)  and  (5),  and  the 
optimal  amplification  factors  were  determined  from  results  for  M2.  However,  theoretically  diurnal 
internal  tides  become  trapped  in  latitudes  higher  than  ~30°  and  no  barotropic  to  baroclinic  energy 
conversion  due  to  freely  propagating  internal  waves  would  result.  Because  we  apply  the  dissipation 
matrix  for  M2,  some  dissipation  of  diurnal  tides  does  occur  at  these  higher  latitudes  in  our  model. 
However,  without  separating  the  modes  and  taking  into  account  the  influences  from  the  other  con¬ 
stituents  similar  to  methods  for  bottom  friction  (e.g  Le  Provost  and  Lyard,  1997),  it  is  unclear 
how  selective  dissipation  for  each  mode  is  possible  using  this  type  of  parameterization  in  a  forward 
model.  Our  assumption  for  internal  tide  dissipation  is  that  the  M2  tidal  wave  dominates  the  signal. 
Nevertheless,  as  shown  in  Fig.  12,  smaller  amplification  factors  are  optimal  for  the  Ki  tidal  wave 
than  those  for  M2.  Furthermore,  in  the  Luzon  Strait  only  resonance  of  the  semi-diurnal  internal 
tides  occur  so  theoretically  the  multiplier  coefficients  Cluzoti  should  not  be  applied  to  the  diurnal 
tides.  Perhaps,  including  information  from  measurements  and  operational  baroclinic  3D  models, 
e.g.  HYCOM  (Chassignet  et  al.,  2007),  may  provide  us  with  an  opportunity  to  locally  improve 
dissipation  matrices  for  depth-integrated  barotropic  models,  although  we  are  still  somewhat  limited 
by  the  assumptions  of  the  underlying  parameterization. 

5.4  Bottom  Friction  Dissipation 

This  section  summarizes  the  effect  of  implementing  the  spatially  varying  Cf  map  (Fig  3)  based  on 
sediment  types  and  tidal  current  speeds  into  the  IndWPac  model.  Firstly,  it  is  useful  to  highlight 
regions  where  we  expect  bottom  friction  to  have  a  large  effect.  Zaron  (2017)  recently  introduced  a 
bottom  friction  number,  Rf  for  the  kth  constituent  to  quantitatively  illustrate  this: 

(cfuf\uk\/h)*  ,20) 

(wfc|ufc|)2  +  (CfUf\uk\/h)2 

where  cuk  is  the  tidal  frequency  of  kth  constituent.  The  term  cjk  corresponds  to  local  acceleration. 
Since  depth  and  velocities  are  highly  correlated  (velocity  and  Cf  are  to  some  extent  also  correlated 
but  much  less  so)  it  is  clear  from  (20)  that  the  effect  of  Cf  is  in  fact  secondary  to  the  effects  of  h . 
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Thus,  a  global  map  of  Rf  based  on  the  spatially  constant  Cf  =  2.5x10  3  simulation  ought  to  suffice 
for  visualization  purposes  (Fig  18). 


Figure  18:  M2  bottom  friction  number  Rf  based  on  the  Comp  +  IT  +  SC  model  setup. 

Rf  shows  a  similar  pattern  for  the  Ki  tidal  wave  (not  shown)  as  the  plotted  M2  tidal  wave.  On 
the  continental  shelves  at  depths  ~100  m,  Rf  is  generally  in  the  0. 1-0.5  range,  and  only  becomes 
larger  than  0.5  close  to  the  coast  in  depths  much  less  than  50  m  (see  also  Zaron,  2017).  Regions 
where  Rf  is  large  correlate  to  areas  where  the  tidal  solutions  will  be  most  impacted  by  any  variability 
in  Cf.  When  implementing  the  spatially  varying  Cf  map  (Fig  3),  specific  areas  with  relatively  large 
Rf  that  are  noticeably  different  to  Cf  =  2.5xl0“3  include:  the  Yellow  Sea  and  southern  Japan 
(small  Cf  in  Yellow  Sea,  large  Cf  just  south  of  Yellow  Sea  and  in  the  Seto  Inland  Sea);  the  area 
between  the  Java  Sea  and  South  China  Sea  near  Singapore  (mostly  large  Cf  with  pockets  of  small 
Cf );  and  Timor  and  Arafura  Seas  (both  large  and  small  Cf).  The  change  in  amplitudes  for  the 
M2  and  Ki  tidal  waves  in  these  three  regions  when  using  the  spatially  varying  Cf  map  over  the 
spatially  constant  Cf  =  2.5xl0-3  are  illustrated  in  Fig  19.  Changes  in  the  RMS  discrepancy  at  the 
coastal  tide  gauges  are  also  plotted.  The  mean  discrepancies  at  the  coastal  tide  gauges  within  the 
boxed  regions  are  summarized  in  Table  5.  The  following  sections  detail  the  findings  at  each  region 
individually,  followed  by  a  discussion  of  the  results. 
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Figure  19:  (a)  M2  and  (b)  Ki  amplitude  differences  due  to  changes  in  bottom  friction  coefficients  Cf  (case 
with  spatially  varying  Cf  map  minus  case  with  spatially  constant  Cf  =  2.5x10  3).  Circles  indicate  the 
change  in  RMS  discrepancies  at  coastal  tide  gauges  (negative  indicates  reduction  in  discrepancy  for  spatially 
varying  Cf  model  setup);  (i)  Yellow  Sea  and  southern  Japan,  (ii)  area  between  the  Java  Sea  and  South 
China  Sea,  (iii)  Timor  and  Arafura  Seas. 
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5.4.1  Yellow  Sea  and  southern  Japan 

In  the  census  sediment  database,  the  Yellow  Sea  (including  Bohai  Sea)  is  designated  as  a  mud 
sediment  type  and  the  tidal  currents  are  very  large  over  the  shallow  basin.  This  leads  to  small 
values  of  Cf  between  l.OxlO"3  and  2.0xl0-3,  except  close  to  the  coastline  where  Cf  becomes  large 
due  to  small  depths  in  (6).  A  patch  of  sand  just  south  of  the  Yellow  Sea  causes  Cf  to  exceed 
4.0xl0-3  here.  Additionally,  a  sand  zone  throughout  most  of  the  Seto  Inland  Sea  induces  a  large 
Cf  (>  4.0xl0-3). 

The  tidal  amplitudes  increase  due  to  the  spatially  varying  Cf  in  most  of  the  Yellow  Sea,  and  de¬ 
crease  just  south  of  the  Yellow  Sea  due  to  the  mud  and  sand  zones  respectively  for  both  constituents. 
The  M2  amplitude  is  also  decreased  in  the  western  Seto  Inland  Sea  but  there  is  no  noticeable  change 
for  Ki.  Impressively,  the  RMS  discrepancy  is  decreased  almost  everywhere  for  both  constituents 
aside  from  a  couple  of  outliers  in  the  Yellow  Sea  and  for  a  group  of  stations  in  the  eastern  Seto 
Inland  Sea  for  M2.  For  example,  the  discrepancy  is  decreased  by  up  to  25  cm  and  5  cm  for  M2  and 
Ki  respectively  at  many  Yellow  Sea  locations.  It  appears  the  the  combination  of  the  small  friction  in 
the  Yellow  and  Bohai  Seas  and  the  large  friction  just  south  of  the  Yellow  Sea  and  in  the  Korea- Japan 
strait  result  in  systematic  positive  changes  to  the  solution.  Overall,  the  mean  RMS  discrepancies 
versus  coastal  tide  gauges  decrease  by  4.9  cm  (28%)  for  M2  and  0.71  cm  (19%)  for  Ki  due  to  the 
spatially  varying  Cf  here. 

5.4.2  Area  between  the  Java  Sea  and  South  China  Sea 

The  region  in  between  the  Java  Sea  and  South  China  Sea  is  predominantly  designated  as  a  sand  zone 
(Cf  >  4.0xl0“3)  with  pockets  of  fine-grain  calcareous  sediment  (Cf  &  2.0xl0“3)  and  volcanic  ash 
(Cf  «  3.0xl0~3)  in  the  census  sediment  database.  There  are  only  a  small  number  of  data-points  in 
the  census  to  back  up  these  sediment  types. 

The  M2  amplitudes  decrease  most  significantly  in  the  region  close  to  Singapore  due  to  the  higher 
values  of  Cf  and  Rf  here.  In  response  to  this  decrease  in  amplitude,  the  M2  RMS  discrepancies 
increase  by  approximately  5  cm  here.  A  large  scale  decrease  in  Ki  amplitudes  occur  southeast  of 
Singapore  but  there  are  no  coastal  tide  gauges  here  to  measure  the  effect  on  the  discrepancy.  A  band 
of  increased  Ki  amplitude  in  between  Singapore  and  Borneo  which  increases  the  RMS  discrepancies 
at  the  tide  gauges  is  also  present.  The  mean  total  free  surface  RMS  discrepancies  are  increased  versus 
tide  gauges  by  1.4  cm  (8.5%)  overall  (Table  5).  Thus,  increasing  Cf  from  the  base  value  of  0.025 
to  the  sediment /current  informed  estimate  clearly  degraded  the  accuracy  in  both  the  semi-diurnal 
and  diurnal  constituents  suggesting  that  either  the  census  sediment  database  is  not  correct  or  that 
the  high  current  speeds  in  the  Malacca  Strait  artificially  increase  the  friction  coefficient  for  the  sand 
sediment  type.  We  believe  that  the  comprehensive  bathymetry  in  this  heavily  trafficked  region  is 
reliable  in  that  it  matches  navigation  charts  well. 

5.4.3  Timor  and  Arafura  Seas 

The  Timor  Sea  is  designated  as  a  sand  sediment  type  in  the  census  sediment  database  while  the 
Arafura  Sea  is  designated  as  mainly  fine-grained  calcareous  sediment  with  a  couple  of  pockets  of 
sand  types.  As  a  result,  in  most  of  the  Timor  Sea  Cf  >  4.0xl0~3  except  for  nearshore  in  the 
Joseph  Bonaparte  Gulf  where  the  tidal  velocities  are  very  large  causing  Cf  to  decrease  here  as  the 
fine-grained  sediment  bedforms  are  washed  out.  The  tidal  velocities  are  fairly  small  in  the  Gulf  of 
Carpentaria  so  even  though  some  of  sediment  is  fine-grained,  Cf  is  similar  to  the  standard  value 
(2.5xl0~3).  North  of  Arnhem  the  sediment  is  fine-grained  and  the  tidal  currents  are  fairly  large 
hence  Cf  becomes  smaller  than  the  standard  value. 
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Because  most  of  the  region  has  larger  values  of  Cf  than  the  standard  value,  tidal  amplitudes 
tend  to  decrease  for  both  constituents  on  the  whole  (Fig  19  (iii)).  Exceptions  are  the  Van  Diemen 
Gulf  (small  Cf ),  north  of  Arnhem  (small  Cf  nearshore  with  high  Cf  just  offshore),  and  east  of 
the  Torres  Strait,  but  only  for  M2.  The  mechanism  for  the  latter  is  mainly  through  dampening  of 
the  energy  fluxes  traveling  east  towards  the  Torres  Strait  as  they  encounter  the  sand  zones  in  the 
Timor  and  Arafura  Seas.  This  results  in  less  resistance  to  the  westward  directed  energy  fluxes  into 
the  Torres  Strait  increasing  amplitudes  to  the  east  of  the  strait  (see  Fig  11).  In  general,  using  the 
spatially  varying  Cf  led  to  poorer  results  for  both  constituents.  Only  deep  in  the  Joseph  Bonaparte 
Gulf  and  at  a  couple  of  stations  near  the  Van  Diemen  Gulf  with  smaller  Cf  values  did  the  M2  RMS 
discrepancy  decrease.  Mean  total  free  surface  RMS  discrepancies  at  tide  gauges  increased  by  1.7  cm 
(6.6%)  due  to  the  spatially  varying  Cf.  Again,  this  could  be  related  to  the  sediment  information 
derived  from  the  census  sediment  data  or  from  the  friction  coefficient  estimation.  In  addition,  the 
degree  of  uncertainty  of  the  bathymetry  is  high  as  explored  in  §5.2,  which  will  substantially  influence 
the  dissipation  and  the  energy  fluxes. 

5.4.4  Discussion 

Results  summarized  in  the  above  sections  highlight  one  region  (the  Yellow  Sea  and  southern  Japan) 
that  experienced  wholesale  decreases  to  the  discrepancies  at  tidal  stations  by  significant  magnitudes, 
and  two  regions  with  small  increases  in  discrepancies,  when  using  the  spatially  varying  Cf  map  over 
the  standard  Cf  =  2.5  xlO-3.  The  Yellow  Sea  region  represents  a  flagship  result  of  the  possibilities 
of  using  a  data-informed  approach  to  estimating  bottom  friction  coefficients  over  using  a  standard 
value.  Furthermore,  the  distribution  of  Cf  used  in  the  region  is  in  close  agreement  to  previous 
studies.  For  example,  Lefevre,  Provost,  and  Lyard  (2000)  found  a  small  value  of  Cf  (1.5xl0-3)  to 
be  optimal  throughout  the  East  China  Sea/ Yellow  Sea  region.  More  interestingly  the  distribution 
of  Cf  closely  resembles  the  optimal  one  determined  through  the  adjoint  method  (Lu  and  Zhang, 
2006).  This  may  appear  to  indicate  that  the  adjoint  method  can  estimate  Cf  values  that  correspond 
to  physical  characteristics  of  the  seabed  such  as  the  muddy  nature  throughout  most  of  the  Yellow 
Sea.  However,  it  is  not  clear  whether  the  sand  zone  south  of  the  Yellow  Sea  creating  an  increase  in 
Cf  similar  to  the  optimal  distribution  of  Cf  in  Lu  and  Zhang  (2006)  was  simply  a  coincidence  or 
not  due  to  the  sparseness  of  data  points  in  the  census  here  and  the  tendency  of  sand  zones  in  the 
other  two  regions  to  in  fact  increase  discrepancies. 

The  Arafura  Sea  in  the  census  is  mainly  fine-grained  calcareous  but  there  are  two  pockets  of 
sand  zones.  However,  according  to  the  grain  size  map  presented  in  Porter-Smith  et  al.  (2004),  we 
expect  the  Arafura  Sea  and  Gulf  of  Carpentaria  to  be  fairly  muddy  with  small  grain  sizes.  In  the 
Timor  Sea  the  census  contains  many  data  points  that  are  used  to  determine  that  the  sediment  is 
sandy  here.  This  may  be  the  case,  but  our  definition  of  the  sediment  grain  size  of  sand  does  not 
seem  to  match  Porter-Smith  et  al.  (2004),  which  indicates  fairly  small  grain  sizes  for  the  region. 
The  rest  of  the  Australian  shelf  has  fairly  large  grain  sizes  (Porter-Smith  et  al.,  2004),  but  they  are 
mainly  designated  as  ooze  types  in  the  census.  It  can  be  assumed  that  there  are  similar  issues  with 
sediment  types  and  the  correlation  of  these  to  physical  sediment  characteristics  in  the  area  between 
the  Java  Sea  and  South  China  Sea.  This  indicates  that  the  census  is  generally  insufficient  for  our 
purpose  due  to  the  relative  sparseness  of  the  data,  particularly  of  terrestrial  type  sediments,  and 
because  it  requires  us  to  estimate  the  grain  size  and  density  from  the  sediment  types.  The  physical 
characteristic  of  sand,  for  example,  can  vary  considerably. 

We  have  experimented  with  changing  some  of  the  sediment  type  definitions  of  the  census.  For 
example,  the  sandy  sediment  definition  in  the  area  between  Java  Sea  and  South  China  Sea  was 
changed  to  fine-grained  calcareous  sediment  which  decreases  the  discrepancies  around  Singapore  for 
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M2.  Similarly,  all  of  the  Arafura  Sea  was  changed  to  a  muddy  sediment  type  and  Timor  Sea  to 
fine-grained  calcareous  sediment  (based  on  the  grain  size  figure  in  Porter-Smith  et  al.  (2004)).  This 
led  to  generally  improved  results.  In  future  studies  if  actual  physical  sediment  characteristics  from 
databases  (e.g  Porter-Smith  et  al.,  2004;  Buczkowski  et  al.,  2006)  could  be  adopted,  with  a  focus  on 
regions  where  Rf  exceeds  0.4-0. 5  in  particular,  there  are  indications  that  non-trivial  improvements 
to  tidal  solutions  can  be  achieved  (e.g.  Yellow  Sea,  Bohai  Sea,  and  Seto  Inland  Sea). 

6  Conclusions 

This  study  presented  a  finite-element  barotropic  model  of  the  Indian  and  Western  Pacific  Oceans 
with  elemental  resolution  ranging  from  as  small  as  100  m  (in  Hong  Kong)  up  to  25  km  in  less 
barotropically  interesting  areas  of  the  deep  ocean.  Most  of  the  resolution  at  the  coast  is  1  km. 
Bathymetry  was  sourced  predominantly  from  the  global  SRTM  bathymetric  database  in  addition 
to  local  high-resolution  datasets  and  hand-edits.  We  began  by  making  comparisons  with  both  the 
data-assimilated  TPX08  atlas  and  tidal  constituents  at  tide  gauges  in  deep,  continental  shelf  and 
slope  waters,  and  at  the  coast.  This  was  followed  by  a  presentation  of  the  sensitivities  to  lateral 
boundary  conditions,  bathymetry,  internal  tide  dissipation,  and  bottom  friction  dissipation.  Within 
each  of  these  sections  a  discussion  of  the  findings  and  implications  with  regards  to  the  sensitivity 
for  that  component  was  presented.  The  key  results  (e.g.  RMS  discrepancies  for  the  Comp  +  IT  + 
SV  model  setup)  and  conclusions  for  each  region  (deep,  shelf,  and  coastal)  follows. 

Deep  water  M2  mean  RMS  discrepancies  against  TPX08  are  2.9  cm  (RMSE  =  3.6  cm)  which 
is  marginally  better  than  reported  for  other  forward  hydrodynamic  models.  However,  the  total  free 
surface  RMS  discrepancies  at  deep  water  tide  gauges  are  2.3  times  those  of  the  TPX08  atlas.  Poorly 
placed  elevation  specified  lateral  boundaries  led  to  global  resonant  amplifications  of  the  lunar  semi¬ 
diurnal  modes.  An  absorption-generation  sponge  zone  suppressed  the  resonant  amplifications  but  it 
relies  on  data-assimilated  model  fluxes  (e.g.  TPX08  and  similar  models)  which  may  not  be  as  reliable 
as  the  elevations  from  these  models.  A  comprehensive  global  forward  model  may  have  advantages 
in  eliminating  the  uncertainties  from  the  boundary  conditions.  Strictly  bathymetric  effects  wTere 
not  of  great  importance  to  the  deep  water  solution,  however  internal  tide  dissipation  (that  relies 
on  topographic  features  and  slopes  for  parameterization)  was  shown  to  be  the  key  control  in  deep 
water.  The  Nycander  method  for  internal  tide  dissipation  was  shown  to  obtain  marginally  superior 
results  to  the  Directional  method,  but  the  latter  was  more  dissipative  in  the  Indian  Ocean  and  hence 
resulted  in  smaller  M2  amplitudes  here  (which  match  TPX08  slightly  better).  There  was  evidence 
that  there  is  scope  to  further  improve  the  deep  water  solution  through  internal  tide  dissipation  but 
it  is  probably  not  possible  using  the  same  paradigm  (strict  reliance  on  the  Nycander  equation  and 
calibration  of  a  global  amplification  factor)  as  presented  here.  Significant  local  modifications  based 
on  3D  baroclinic  models  and  measurements,  including  improvements  to  internal  tide  dissipation  in 
shallower  waters  are  likely  required. 

Continental  shelf  and  slope  water  M2  mean  RMS  discrepancies  against  TPX08  are  6.5  cm  (RMSE 
=  10.1  cm).  This  was  shown  to  be  significantly  superior  to  those  reported  for  other  forward  hy¬ 
drodynamic  models  (RMSE  =  19-28  cm  (Stammer  et  al.,  2014),  albeit  these  are  global  errors,  the 
total  energy  density  is  at  least  as  large  in  the  IndWPac  domain  compared  to  the  rest  of  the  world’s 
ocean).  One  of  the  most  important  factors  for  improvement  was  shown  to  be  the  inclusion  of  local 
high-resolution  bathymetry.  Notable  changes  in  the  M2  amplitude  and  a  corresponding  reduction 
in  the  RMS  discrepancies  against  TPX08  and  at  tide  gauges  were  evident  in  the  greater  Yellow  Sea 
region  due  to  changes  in  bottom  friction  coefficients  based  on  combinations  of  muddy  and  sandy 
sediment  types.  However,  the  regions  around  the  Java  and  South  China  Seas,  and  Timor  and  Ara- 
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fura  Seas  did  not  generate  significant  nor  positive  changes  to  the  discrepancy.  If  more  complete 
databases  of  physical  characteristics  of  sediment  are  made  available,  combined  with  accurate  local 
bathymetric  data  it  may  be  possible  to  improve  solutions  elsewhere,  particularly  in  resonant  basins 
(e.g.  Gulfs  of  Khambhat  and  Kutch  on  the  west  coast  of  India;  King  Sound,  Joseph  Bonaparte  Gulf 
and  Van  Diemen  Gulf  in  northern  Australia,  among  others)  due  to  spatially  varying  bottom  friction 
coefficients.  Nevertheless,  we  are  still  limited  by  residual  discrepancies  from  deeper  waters  and  the 
uncertainties  of  internal  tide  dissipation  in  shallower  waters. 

The  discrepancies  at  the  coast  for  the  IndWPac  model  were  not  significantly  different  from  those 
obtained  further  offshore  on  the  shelf.  Furthermore,  the  mean  total  free  surface  RMS  discrepancies 
at  coastal  tide  gauges  ( Dtg  =  14  cm)  were  2.6  cm  smaller  than  those  of  the  TPOX8  atlas.  However, 
the  discrepancy  at  the  majority  of  locations  are  smaller  for  the  TPX08  atlas  due  to  data-assimilation 
offshore  and  at  selected  tide  gauges.  The  large  mean  RMS  discrepancy  for  the  TPX08  atlas  is  likely 
related  to  TPX085s  coarser  resolution  not  resolving  certain  nearshore  features  and  harbor  complexes 
in  detail  and  therefore  not  correctly  propagating  the  tides  into  them.  For  example,  bathymetric  and 
bottom  friction  controls  were  both  found  to  play  a  very  important  role  nearshore,  and  in  some  cases 
dominate  the  reasons  for  discrepancy  due  to  resonance  in  a  basin  or  inadequate  connectivity  into  a 
bay.  In  contrast,  the  IndWPac  model  does  not  have  such  a  large  number  or  magnitude  of  outlier 
locations.  These  results  are  an  indication  that  the  model  adequately  captures  a  large  amount  of 
the  nearshore  physics  throughout  the  domain.  Thus,  the  model  is  potentially  suitable  to  simulate  a 
great  range  of  shallow  water  physics  within  the  region,  specifically  into  detailed  harbor  complexes 
and  other  nearshore  features  where  the  tide  gauges  are  located. 

If  we  are  solely  interested  in  tidal  elevations,  then  the  simple  answer  is  to  use  data-assimilation 
within  the  IndWPac  model  to  achieve  highly  accurate  solutions,  from  the  deep  ocean  all  the  way  to 
the  well-resolved  coastal  regions.  However,  in  many  other  applications,  such  as  the  forecasting  and 
analysis  of  coupled  surge,  tide  and  wave  processes,  capturing  the  large-scale  responses  to  meteorology, 
and  modeling  the  shallow  water  physics  including  the  nonlinear  interactions  of  the  processes  becomes 
vital.  In  order  to  accomplish  this,  correctly  specifying  high-resolution  bathymetry  and  topography 
become  controlling  factors.  Furthermore,  physics  based  improvements  to  more  accurately  quantify 
dissipation  within  forward  barotropic  models  are  possible  offshore,  through  coupling  to  coarser  3D 
baroclinic  numerical  models,  and  nearshore,  through  bottom  bedform  and  sediment  roughness  data. 
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A  Effective  Sediment  Roughness  Equations 


The  equations  for  the  effective  sediment  roughness  ks  are  taken  from  Rijn  (2007).  ks  is  calculated 
from  the  vector  sum  of  roughnesses  from  the  different  bedform  types: 

ks  =  (k2sr  +  k2m  +  k2d )  (21) 

where  kS)r ,  /cSjm,  ks,d  are  the  ripple,  mega-ripple,  and  dune  related  roughnesses  respectively.  Equa¬ 
tions  for  each  rely  on  the  current  mobility  parameter  0: 


if  — 
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2 
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(s  -  l)gdb0 


(22) 


where  uj  is  the  effective  mean  current  speed  (7),  s  =  ps/po  is  the  relative  sediment  density  ( ps  is 
the  sediment  density),  g  is  the  acceleration  due  to  gravity,  and  d50  is  the  median  sediment  grain 
diameter.  The  equations  for  each  individual  roughness  type  are  then: 


ks,r  =  fscdso  (85  -  65  tanh[O.O15(0  -  150)])  (23) 

ks,m  “  max  (min  (0.02, 200c?5o)  >  2 e~5ffsh  [1  —  exp(— 0.05*0)]  (550  —  0))  (24) 

kSid  =  max  (0, 8e~5ffsh  [1  -  exp(-O.O20)]  (600  -  0))  (25) 
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where  h  is  the  still  water  depth,  fsc  is  the  “factor  which  expresses  the  effect  of  a  gradually  decreasing 
ripple  roughness  for  very  coarse  sediment  beds” ,  and  ffs  is  the  “factor  which  expresses  the  effect  of 
a  gradually  decreasing  mega-ripple  roughness  for  very  fine  sediment  beds”  (Rijn,  2007): 


fca  ~  min 
ffa  =  min 


f  0.25 djgrav  \ 

’  V  *0  ) 
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1.5  dsan(i 


(26) 

(27) 


where  dgrav  =2x10  3  m,  and  dsand  =  6.2x10  5  m.  Here,  all  sediments  are  assumed  to  have 
g?5q  >  dsiit  =  3.2xl0“5  m,  in  which  otherwise  a  lower  limit  of  ks  =  20 dsut  is  applied. 
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